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ABSTRACT 


This  paper  describes  the  development  of  a  6-degree  of  freedom  (6-DOF), 
nonlinear,  miniature  rotary-wing  unmanned  aerial  vehicle  (RW  UAV)  simulation 
environment  using  MathWorks  Simulink  simulation  software.  In  addition  to  the  modeling 
process,  this  research  also  conducts  flight-path  controller  design  using  Proportional- 
Derivative  (PD)  control  techniques.  This  model’s  development  is  motivated  by  the  desire 
to  enable  a  rapid  prototyping  platform  for  design  and  implementation  of  various  flight 
control  techniques  with  further  seamless  transition  to  the  hardware  in  the  loop  (HIL)  and 
flight-testing.  The  T-Rex  Align  600  remote  controlled  helicopter  with  COTS  autopilot 
was  chosen  as  a  prototype  rotary  UAV  platform. 

The  development  of  the  nonlinear  simulation  model  is  implemented  starting  with 
extensive  literature  review  of  helicopter  aerodynamics  and  flight  dynamics  theory  and 
applying  the  mathematical  models  of  the  helicopter  components  to  generate  helicopter 
inertial  frame  motion  simulations  from  operator  commands.  The  primary  helicopter 
components  modeled  in  this  thesis  include  the  helicopter  main  rotor  inflow,  thrust, 
flapping  dynamics,  as  well  as  the  tail  rotor  inflow  and  thrust  responses.  The  inertial  frame 
motions  are  animated  using  the  Flight  Gear  Version  0.9.8  software. 

After  obtaining  simulations  with  verifiable  results,  the  nonlinear  model  is 
linearized  about  the  hovering  flight  condition  and  a  linear  model  is  extracted.  Lastly,  the 
PD  controller  is  designed  and  flight  path  software  in  the  loop  (SIL)  test  results  are 
presented  and  explained.  The  SIL  tests  are  conducted  for  autonomous  flight  along 
specified  rectangular  and  figure-8  flight  paths. 
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EXECUTIVE  SUMMARY 


The  use  of  rotary-wing  unmanned  aerial  vehieles  (RW  UAV)  brings  an 
unparalleled  advantage  in  defense  applieations.  Unlike  its  fixed-wing  eounterparts, 
rotary- wing  platforms  provide  mueh  more  flexibility  and  maneuverability  in  the  joint 
battle  spaee.  The  RW  UAV  is  more  flexible  beeause  it  does  not  require  a  eleared  runway 
or  additional  airspaee  for  landing  or  takeoff  maneuvers.  Beeause  of  its  vertieal  takeoff, 
landing,  and  flight  eapabilities,  the  RW  UAV  makes  an  exeellent  platform  for  low 
altitude  aerial  reeonnaissanee  of  urban  terrain. 

Harnessing  this  potential  demands  the  development  of  autopilot  eontrol 
algorithms  that  are  testable  and  robust.  Autopilot  and  RW  UAV  hardware  is  expensive, 
and  no  eurrent  testing  software  ean  simulate  the  flight  dynamies  of  the  Meehanieal  and 
Astronautieal  Engineering  (MAE)  Department’s  miniature  RW  platform,  the  T-Rex 
Align  600.  The  development  of  a  simulation  environment  is  a  step  eloser  towards  being 
able  to  design  and  test  eontrol  algorithms  while  lowering  the  risk  of  breaking  expensive 
hardware  eomponents  during  flight  tests.  The  intent  of  the  thesis  research  conducted 
herein  is  to  enable  and  facilitate  the  development  of  a  miniature  RW  simulation 
environment.  It  is  hoped  that  the  results  provided  here  will  spawn  further  interest  and 
development  efforts  that  will  lead  to  a  more  accurate  and  reliable  control  test  bed. 

The  first  stage  in  the  development  of  the  simulation  environment  is  the  modeling 
of  the  6-DOE  RW  dynamics  of  the  miniature  helicopter.  This  is  done  by  conducting 
extensive  literature  review  of  the  authoritative  helicopter  aerodynamics  and  flight 
dynamics  theory,  and  applying  the  mathematical  models  of  the  helicopter  components  to 
generate  helicopter  inertial  frame  motion  simulations  from  control  inputs.  The  primary 
helicopter  components  modeled  in  this  thesis  include  the  helicopter  main  rotor  inflow, 
thrust,  flapping  dynamics,  as  well  as  the  tail  rotor  inflow  and  thrust  responses.  The  RW 
model  responses  are  then  qualitatively  compared  to  the  data  and  results  of  [1]  through 
[7].  The  software  6-DOE  model  was  created  using  MathWorks  Simulink  simulation 
software.  The  inertial  frame  motions  are  animated  using  the  Elight  Gear  Version  0.9.8 
software. 


The  second  stage  in  the  development  of  the  simulation  environment  is  the 
linearization  of  the  nonlinear  model.  Controller  design  focuses  on  the  hovering  flight 
regime  and  linearization  helps  to  uncover  the  state  boundaries  for  which  the  controller 
will  be  useful.  The  last  stage  of  the  development  of  the  simulation  environment  was  the 
proportional-derivative  (PD)  controller  tuning  and  implementation  on  the  nonlinear 
model.  The  PD  controller  performance  was  tested  by  analyzing  flight  path  tracking 
performance. 

The  complete  6-DOF  model  performs  well  in  modeling  the  dynamics  of  a 
miniature  RW  helicopter  as  compared  to  the  predictions  and  flight  test  data  provided  by 
[1]  through  [7].  In  addition,  analysis  of  the  flight  path  tracking  performance  leads  to  the 
conclusion  that  the  derived  PD  gains  provide  satisfactory  control  of  the  model.  Maximum 
errors  in  position  and  yaw  angle  response  are  expected  and  are  within  tolerance. 

It  is  hoped  that  the  research  and  results  from  this  effort  will  lead  to  further 
development  and  refinement  of  the  simulation  environment  that  has  been  developed  up  to 
this  point. 
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NOMENCLATURE 


Latin  Variables 


Cj  Coefficient  of  thrust  [•] 

T 

C  = - 

^  p{QRf  ttR' 

4z  Helicopter  Principal  Moments  of  Inertia  [.18  .34  .28]  [kg  m  ] 

Kp  MR  Hub  Stiffness  [Nm/rad] 

L,  M,  N  Helicopter  Roll,  Pitch,  and  Yaw  moments  [N  m] 

Lb,  Ma  Roll  Moment  and  Pitch  Moment  Derivative  [Nm/rad] 

Tail  rotor  moment  arm.  Vertical  length  from  tail  rotor  hub  to  .08[m] 

helicopter  c.g.  line  of  sight 

Tail  rotor  moment  arm.  Length  from  helicopter  c.g.  to  tail  rotor  .91  [m] 

hub 

• 

m  Mass  flow  rate  of  air  through  rotor  [kg/sec] 

N  Number  of  blades 

p,  q,  r  Helicopter  body  roll,  pitch,  and  yaw  rates,  respectively  [rad/sec] 

T  Main  rotor  thrust  [N] 

Ttr  Tail  rotor  thrust  [N] 

u,  V,  w  Helicopter  airspeed  in  the  x,  y,  and  z  directions,  respectively  [m/s] 

Utr,  Vtr,  Wtr  Tail  rotor  hub  airspeed  in  the  x,  y,  and  z  directions,  respectively  [m/s] 

V-  Normal  component  of  rotor  induced  velocity  [m/s] 

y  „  Induced  flow  of  air  through  rotor,  far  downstream  of  rotor  [m/s] 

V,  V„  Helicopter  airspeed  (magnitude)  [m/s] 

Vn  Velocity  of  air  normal  to  rotor  [m/s] 

X,  Y,  Z  Helicopter  forces  in  the  X,  Y  and  Z  directions  [N] 
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Greek  Variables 


a  Angle  of  attack.  Angle  between  horizon  and  main  rotor  tip  path  [rad] 

Plane 

Plan  ’  Plat  Longitudinal  and  Lateral  Main  Rotor  Flap  Angles  [rad] 

3  3 

p  Air  density  (1.275  kg/m  at  sea  level)  [kg/m] 

/I,  ,  Inflow  ratio  for  main  rotor  [•] 

p  Advance  ratio.  Ratio  between  helicopter  translational  velocity  [•] 

and  main  rotor  tip  speed 

Advance  ratio  along  the  z  axis.  Ratio  of  vertical  velocity  to  [•] 

main  rotor  tip  speed 

Collective  lever  input  ( 0  [•] 

Lateral  stick  input  (left-right  -  [*] 

Longitudinal  stick  input  (forward-aft  -  1  ^  <  1 )  [•] 

Rudder  pedal  input  (left-right  -  [•] 

Q  Main  Rotor  rotation  rate  167[rad/sec] 

Tail  Rotor  rotation  rate  778.2  [rad/sec] 

6^  Collective  pitch-angle  of  main  rotor  blade  [rad] 

Lateral  cyclic-pitch  angle 

Angle  of  main  rotor  blade  when  blade  is  aligned  with  the  +y  -axis  [rad] 
Longitudinal  cyclic-pitch  angle.  Angle  of  main  rotor  blade  when 
blade  is  aligned  with  the  -x  -axis  [rad] 

Collective  pitch-angle  of  tail  rotor  blade  [rad] 

(j),  6,  y/  Euler  orientation  angles  [rad] 

z Main  Rotor  Time  Constant  [sec] 
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Constants 


a 

Main  rotor  lift  curve  slope 

5.5[rad'^] 

atr 

Tail  rotor  lift  curve  slope 

5.0[rad-^] 

c 

Rotor  blade  chord 

[m] 

Ad 

Main  rotor  disk  area 

1.887[m^] 

g 

Gravitational  acceleration 

9.82[in/s^] 

m 

Mass  of  helicopter 

8.2[kg] 

R 

Main  Rotor  radius 

.775  [m] 

Rtr 

Tail  Rotor  radius 

.13[m] 

s 

Rotor  Solidity  Nc/(7iR) 

[•] 

Abbreviations 

MR  Main  Rotor 

fus  Fuselage 

tr  Tail  Rotor 
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I.  INTRODUCTION 


A.  MOTIVATION 

Rotary- wing  Unmanned  Aerial  Vehicles  (RW  UAV)  provide  capabilities  that 
fixed-wing  UAVs  cannot.  Unlike  fixed-wing  UAVs,  RW  UAVs  provide  enhanced 
flexibility  in  the  joint  Navy,  Marine  Corps,  Army,  and  Air  Force  battle  space  by  enabling 
flexibility  in  aerial  reconnaissance.  In  contrast  to  the  fixed-wing  UAVs,  RW  UAVs 
provide  vertical  take-off  and  landing  (VTOL),  thereby  foregoing  the  necessity  for  cleared 
terrain.  Additionally,  RW  UAVs  provide  the  ability  to  conduct  personnel  clearing  in 
urban  terrain  to  a  degree  that  is  not  possible  with  fixed-wing  UAVs. 

The  advantages  provided  by  a  RW  UAV  make  it  highly  desirable  for  DoD 
applications.  Development  of  a  nonlinear  RW  UAV  software  model  will  facilitate  future 
development  of  RW  UAV  control  design  at  NFS. 

B,  RESEARCH  OBJECTIVES 

The  overall  objective  of  the  research  conducted  for  this  thesis  is  to  provide  a 
simulation  environment  that  will  facilitate  and  enable  future  work  such  as  Hardware  in 
the  Loop  (HIL)  testing  of  the  Micro  Pilot  autopilot  (MP  AP)  and  flight  testing  of  the  MP 
AP  on  the  physical  plant.  Micro  Pilot  software  exists  for  testing  of  controller  code; 
however,  this  software  is  not  specific  to  the  MAE  Department’s  UAV  model  (T-Rex 
Align  600),  nor  is  it  configurable  to  any  other  miniature  RW  platform.  The  simulation 
environment  will  simulate  the  6-Degree  of  Freedom  (6-DOF)  model  of  the  MAE 
Department’s  RW  UAV  and  will  provide  a  test  environment  that  is  configurable  to  any 
specified  model. 

This  research  effort  focuses  on  the  implementation  of  existing  helicopter 
dynamics  theory  to  develop  a  non-linear  6-DOF  RW  UAV  model  using  the  MathWorks 
Simulink  modeling  environment.  Emphasis  is  placed  on  obtaining  results  that  are  in 
agreement  with  the  available  literature  covering  helicopter  dynamics  theory. 
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c. 


T-REX  ALIGN  600  MODEL  DESCRIPTION 


The  T-Rex  Align  600  is  a  remote  eontrolled  model  helieopter.  The  main  rotor  is 
eomposed  of  two  rotor  blades,  made  of  earbon  fiber,  and  a  stabilizer  bar  for  aetive  rate 
damping.  The  T-Rex  Align  uses  an  eleetrie  power  system  and  a  lightweight  carbon  fiber 
frame.  The  canopy  is  made  of  fiberglass  and  the  tail  rotors,  like  the  main  rotors,  are  made 
of  carbon  fiber  material.  Figure  1  illustrates  the  side  view  of  the  T-Rex  Align  600  air 
frame. 


Figure  1 .  T-Rex  Align  600  Model  (Side  View) 
The  T-Rex  Align  600  plant  properties  are  outlined  in  Table  1 . 


PROPERTY 

QUANTITY 

Flying  Mass/  Weight 

3  kg  (29.4  N) 

Main  Rotor  Radius 

.6  m 

Main  Rotor  Diameter 

1.35  m 

Tail  Rotor  Diameter 

.24  m 

Table  1 .  T-Rex  Align  600  Plant  Properties 
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D,  APPROACH 


Helicopter  dynamics  is  dominated  by  the  forces  and  moments  generated  by  the 
helicopter  main  rotor  and  tail  rotor.  First,  pilot  input  to  rotor  blade  responses  are  analyzed 
and  modeled.  Second,  rotor  aerodynamics  is  considered  in  order  to  develop  predictable 
rotor  thrust  responses.  Also,  main  Rotor  flapping  dynamics  are  analyzed  and  modeled. 
The  main  rotor  flapping  is  what  generates  translational  acceleration  and  angular  rotations 
of  the  helicopter  body.  Once  the  main  rotor  thrust  and  flapping  is  modeled,  all  forces  and 
moments  are  summed  and  applied  at  the  helicopter  body  center  of  gravity,  and  vehicle 
motion  is  simulated  using  MathWorks  Simulink  simulation  environment.  Flight  Gear 
software  is  incorporated  with  the  Simulink  model  to  provide  motion  animation. 

There  is  extensive  documentation  of  previous  modeling  efforts  that  have  been 
conducted.  The  approach  of  this  thesis  effort  relies  on  the  modeling  of  known  parameters 
for  the  MIT  instrumented  X-Cell  SE  helicopter  used  in  [1]  and  [2].  The  intent  is  to  model 
a  helicopter  similar  to  the  ME  Department’s  T-Rex  Align  in  order  to  provide  a  platform 
from  which  future  work  in  system  identification  can  uncover  the  specific  model 
parameters  for  this  particular  UAV.  Model  parameters  for  the  MIT  X-Cell  SE  are  listed 
in  Appendix  A. 

After  completion  of  the  modeling,  the  nonlinear  model  is  linearized  about  the 
hovering  flight  condition  in  order  to  uncover  the  state  boundaries  about  which  the 
controller  will  be  useful.  After  linearization,  the  PD  gains  are  estimated,  tuned,  and  tested 
for  flight-path  tracking  performance.  The  flight  paths  that  the  model  will  be  expected  to 
fly  are  the  square  and  figure-8  flight  paths.  Einally,  the  flight  path  tracking  performance 
results  are  discussed. 
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II.  MODELING  INTRODUCTION 


A.  BODY-FIXED  FRAME  AND  HELICOPTER  VARIABLES 


The  modeling  of  the  RW  UAV  uses  only  one  body  eentered  referenee  frame. 
Figure  2  illustrates  the  helieopter  body-fixed  referenee  frame  with  the  origin  at  the  eenter 
of  gravity.  The  variables  are  shown  on  the  body  axes  x,  y,  z .  The  x  -axis  is  positive 
pointing  away  from  the  nose  of  the  UAV,  the  y  -axis  is  positive  towards  the  starboard 
side,  and  the  z  -axis  is  positive  downward.  The  variables  inelude  the  body  veloeities 
u,  V,  w ,  the  Euler  angles  (p,  6,  ,  and  the  body  rates  p,  q,  r . 


hub  plane 

y  _  _  ^ 


.^1 


y.  V,  0.  q 


c.g.  '-■ 


a:,  u,  (j),  p 


I  z,  w,  \u.  r 

t 


Figure  2.  Helieopter  Body-Referenee  Frame  (After  [2]) 

Positive  rotations  are  defined  using  the  right-hand  rule,  with  the  thumb  pointing  in 
the  positive  axis  direetion.  The  rotation  is  defined  by  eurling  the  other  four  fingers  in  the 
direetion  of  rotation.  For  example,  with  the  right  thumb  pointing  in  the  positive  z 
direetion  (pointing  down),  a  positive  yaw  angle  is  defined  by  a  eloekwise  body  rotation. 

The  main  rotor  disk  ean  tilt  about  the  rotor  hub  along  two  degrees  of  freedom, 
about  the  lateral  (y-axis)  and  longitudinal  (x-axis)  direetions.  The  flap  angles  and 

represent  the  main  rotor  flap  angles  about  the  lateral  and  longitudinal  direetions. 
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respectively.  A  positive  lateral  flap  angle  occurs  when  the  main  rotor  tilts  towards 
the  positive  y  -axis  because  this  corresponds  to  a  rotation  in  the  direction  of  positive  roll 
angle  (^).  A  negative  longitudinal  flap  angle  occurs  when  the  main  rotor  tilts 

towards  the  positive  x-axis  because  this  corresponds  to  a  negative  pitch  angle  (6*).  Note 
that  a  positive  pitch  angle  ( ^)  is  defined  by  placing  the  right  hand  face  up,  with  the  right 
hand  thumb  pointing  towards  the  right  (+y  direction).  Curling  the  fingers  in  towards  the 
hands  defines  the  positive  pitch  rotation. 

B,  RIGID-BODY  EQUATIONS  OF  MOTION 

The  helicopter  is  able  to  simultaneously  translate  and  rotate  about  all  6-DOF.  The 
equations  of  motion  for  the  fuselage  six  degrees  of  freedom  are  derived  by  Gavrilets  in 
[1].  Summing  the  moments  and  forces,  due  to  the  main  rotor,  tail  rotor,  and  aerodynamic 
forces  on  the  fuselage,  gives  the  three  translational  equations  of  motion  as  functions  of 
the  UAV  velocity  states  (w,  v,  w)  and  the  forces  acting  on  the  body  center  of  gravity 

i^MR’  ^fus’  for  example)  [1,  p.  33]; 


(^MR+^fus+^,r)  .  , 

u  =  vr-wq-\- - - - gsmO 


m 


v  =  wp-ur-\ - - - +  g  cos  6'sin  (p 


m 


(2.1) 

(2.2) 


i^MR  ^fus  ^tr)  „ 

w-uq-vp  + - - - +  g  cos  &  cos  (, 


m 


(2.3) 


The  X  term,  for  example,  denotes  a  force  in  the  x  direction,  while  the 
subscripts  MR ,  fus  and  tr  represent  the  forces  in  the  x  direction  due  to  the  main  rotor, 
fuselage  and  tail  rotor,  respectively.  The  term  g  is  the  gravity  force. 
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Assuming  the  cross  products  of  inertia  are  small  gives  the  three  equations  of 
motion  describing  the  body-angular  motion  as  a  function  of  the  body  rotation  rates 
{p,  q,  r)  and  moments,  acting  on  the  body  center  of  gravity; 


•  _  qrjlyy  -  ^  (l^r  +  ) 

I  1 

XX  XX 

•  _  -  fxc)  I  ) 

^  I  I 

yy  yy 

•  P^l^Ixx-Iyy) 

r- - - L 

^  zz  ^  zz 

All  of  the  translational  accelerations,  velocities,  angular  rates,  and  Euler  angles 
will  be  generated  using  the  Simulink  6-DOF  block  with  the  main  rotor,  tail  rotor,  drag, 
and  aerodynamic  forces  and  moments  as  inputs.  Figure  3  illustrates  the  forces  and 

moments  acting  on  the  UAV  body.  Note  that  the  main  rotor  thrust 

and  the  drag  force  on  the  body  isF^,„  =  ^fus\  ■  The  tail  rotor  thrust  is 

modeled  to  act  normal  to  the  fin  plane. 


(2.4) 

(2.5) 

(2.6) 


Figure  3.  Helicopter  Forces  and  Moments  (After  [2]) 
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The  outputs  to  the  6-DOF  block  will  be  the  center  of  gravity’s  position,  velocity, 
Euler  angles,  and  Euler  rates.  These  are  derived  by  the  6-DOE  block  by  implementing  the 
state  dynamic  equations  (2.1)  through  (2.6). 

In  summary,  the  model  has  14  vehicle  states; 

x  =  [x  y  z  ^  e  y/  u  V  w  p  q  r 

with  a  body-fixed  coordinate  system  at  the  center  of  gravity.  Positive  rotations  are 
defined  using  the  right  hand  rule  and  the  forces  and  moments  are  summed  at  the  vehicle 
center  of  gravity.  The  Simulink  6-DOE  block  is  used  to  generate  the  vehicle  state 
dynamics. 

The  following  chapter  will  discuss  the  modeling  of  the  necessary  components  that 
will  generate  the  forces  and  moments  on  the  vehicle  center  of  gravity,  as  well  as  the 
modeling  of  the  channel  for  pilot  input  to  control  surface  command. 
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III.  ROTARY-WING  UNMANNED  AERIAL  VEHICLE  (RW  UAV) 

MODELING 

A,  CONTROL  INPUT  TO  CONTROL  SURFACE  COMMAND  DESIGN 

There  are  four  inputs  that  a  pilot/operator  can  give  to  the  UAV.  The  first  input  is 
collective  and,  in  full-scale  helicopters,  is  provided  through  a  lever.  This  lever  is  an 
elongated  version  of  what  looks  like  the  emergency  brake  in  a  car.  A  change  in  the  lever 
position  will  provide  a  change  in  the  main  rotor  collective,  or  the  main  rotor  blade  angle 
of  attack.  This  is  accomplished  through  a  variable  orientation  swash  plate.  As  the  pilot 
provides  collective  input,  the  swash  plate  changes  accordingly.  A  collective  lever  input 
causes  the  swash  plate  to  rise  up  or  down  while  maintaining  a  flat  orientation.  In  the 
UAV  controller,  collective  is  controlled  through  one  of  the  channel  stick  inputs,  normally 
set  at  full  down  by  an  internal  spring  force. 

The  second  type  of  input  available  to  the  pilot/operator  is  the  cyclic.  In  a  full- 
scale  helicopter,  this  is  provided  by  the  stick  immediately  to  the  front  of  the  pilot,  which 
can  rotate  by  varying  degrees.  The  cyclic  input  can  be  classified  by  longitudinal  cyclic 
and  lateral  cyclic.  Longitudinal  cyclic  is  a  cyclic  stick  input  that  affects  the  translational 
acceleration  of  the  vehicle  along  the  x  -axis.  Longitudinal  cyclic  commands  also  affect 
the  pitch  orientation  of  a  helicopter.  This  is  achieved  by  providing  a  forward  or  aft  stick 
command.  A  full  forward  stick  command  causes  the  main  rotor  blade  angle  to  be  a 
maximum  when  it  reaches  the  -x  -axis  and  cycles  through  to  a  minimum  angle  of  attack 
when  the  blade  reaches  the  +x  -axis.  In  a  full-scale  helicopter,  this  is  achieved  by  a  tilting 
of  the  swash  plate,  which  in  turn,  causes  cyclic  (variable)  blade  angle  commands  as  the 
rotor  blade  rotates  around  its  circular  trajectory. 

Lateral  cyclic  is  a  cyclic  stick  input  that  affects  the  translational  acceleration  of 
the  vehicle  along  the  y  -axis.  Lateral  cyclic  also  affects  the  helicopter  roll  orientation. 
This  is  achieved  by  providing  left  and  right  stick  commands.  A  full  right  stick  command 
causes  the  main  rotor  blade  angle  to  be  a  maximum  when  the  rotor  blade  reaches  the  -y  - 
axis;  as  the  blade  cycles  to  the  -l-y-axis,  along  its  circular  trajectory,  the  blade  angle 
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becomes  increasingly  smaller  until  it  reaches  a  minimum  at  the  +y  -axis.  The  maximum 
blade  angle  at  the  -y  position  is  what  generates  a  higher  lift  force  on  the  blade  at  this 
position,  when  compared  to  the  lift  of  the  blade  at  the  +y  position.  This  difference  in  lift 
is  what  causes  the  main  rotor  disk  plane  to  flap.  As  the  rotor  blade  achieves  its  maximum 
angle  of  attack  at  the  -y  position,  the  main  rotor  tilts  in  the  +y  direction,  and  this  causes 
the  helicopter  to  accelerate  in  the  +y  direction.  Again,  this  is  physically  achieved  through 
the  mechanical  tilting  of  the  swash  plate,  which  in  turn  controls  the  cyclic  behavior  of  the 
rotor  angle  of  attack. 

The  last  form  of  pilot  command  is  the  rudder  pedal.  The  rudder  pedal  provides 
collective  input  commands  to  the  tail  rotor,  much  like  the  collective  lever  provides 
collective  commands  to  the  main  rotor.  The  thrust  generated  by  the  tail  rotor  is  necessary 
for  canceling  the  torque  on  the  inertial  frame  that  is  generated  by  the  main  rotor. 

Table  2  outlines  the  four  different  command  types  and  their  physical  meaning  on 
a  handheld  R/C  controller. 


Command 

(Simulink  Pilot  Joystick  Blockset  ID) 

Stick  Input  Value 
Range 

Physical  Meaning  on  Handheld 
Controller 

Collective  Input  (Throttle) 

[0,1] 

[Full  throttle  down,  FT  up] 

Longitudinal  Cyclic  (Pitch) 

[-1,1] 

[Full  Back  Stick,  Full  Fwd  Stick] 

Lateral  Cyclic  (Roll) 

[-1,1] 

[Full  Left  Stick,  Full  Right  Stick] 

Rudder  Collective  (Y aw) 

[-1,1] 

[Full  Left  Pedal,  Full  Right  Pedal] 

Table  2.  Joystick  Block  Value  Range  Outputs 


The  full  positions  provide  maximum  acceleration  that  corresponds  to  the  specified 
input.  Prouty  [3]  outlines  several  helicopters  and  their  collective,  cyclic,  and  tail  rotor 
angle  of  attack  ranges.  These  are  listed  in  Table  3  and  provide  references  for  choosing 
initial  estimates  for  the  model  collective  and  cyclic  ranges. 
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Aircraft 

MR  Collective 

TR  Collective 

Lateral  Cyclic 

Longitudinal  Cyclic 

Agusta  A 109 

o 

0 

1 

ON 

0 

-7°/21° 

±6.25° 

-10.5°/12.5° 

Bell  206L-3 

6° -22° 

-12.45°/19.5° 

±10° 

±10° 

Bell  412 

0 

NO 

1 

o 

O 

-10.15°/19.85° 

-12.6°/6.7° 

±11° 

MBB  BO  105  LS 

-.2°/15° 

-8°/ 20° 

-5.7°/ 4.2° 

-ll°/6° 

Robinson  R22  Beta 

1.5°-14.5° 

-10.6°/19.5° 

-9.5°/ 6° 

-9°/ 11° 

Table  3.  Sample  of  Full-Scale  Aircraft  Control  Ranges  (After  [3]) 


The  response  ranges  that  are  used  for  the  T-Rex  Align  model  are  listed  below; 

X-Cell  60  Simulink  Model  Rotor  Response  Ranges 
MR  Collective  TR  Collective  Lat  Cyclic  Long  Cyclic 

-5°/ +15°  -20°/ +25°°  ±20°  ±20° 

The  Simulink  model  is  designed  so  that  no  input  is  required  to  maintain  the 
vehicle  in  hovering  trim.  Also,  the  cyclic  limits  are  higher  than  the  average  full  scale 
limit  since  the  miniature  UAV  is  more  maneuverable  and  aerobatic. 

The  first  aspect  of  the  model  to  be  designed  is  the  transformation  that  takes  place 
from  the  stick  control  input  [-1,1],  for  example,  to  the  control  surface  command.  The 
stick  input  with  signal  -1  corresponds  to  a  full  left  stick  input  while  a  control  input  signal 
of  1  corresponds  to  a  full  right  stick  input;  a  control  signal  of  0  corresponds  to  the  neutral 
position. 

The  stick  position  [-1,1]  is  coded  in  the  delay  time  between  pulses  of  the  remote 
controller’s  baseband  signal.  The  radio  signal  illustrated  in  Figure  4  corresponds  to  the 
delay  time  for  each  of  five  stick  positions.  The  long  delay  at  the  end  of  the  signal  is  a 
reset  for  the  next  set  of  pulses.  The  receiver  onboard  the  UAV  splits  the  information  from 
each  input  and  sends  it  to  the  corresponding  servo,  in  turn.  Each  of  the  servos  receive  as 
inputs  a  pulse  width  modulated  signal,  with  duty  cycle  information,  typically  in  the  range 
from  600  /us  to  2400  /us .  The  servo  then  decodes  the  duty  cycle  information  and 
generates  a  mechanical  rotation  from  -45° to  +45°.  The  servo’s  mechanical  output  is 
proportional  to  the  duty  cycle.  A  duty  cycle  input  to  the  servo  of  600  /us ,  for  example, 
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will  give  a  servo  output  of  -45° ;  a  duty  cycle  input  of  1,520  jus  will  give  an  output  of 
0°,  and  a  duty  cycle  input  of  2,400  jis  will  give  an  output  of  +45°.  All  angles  in 
between  can  be  achieved  by  varying  the  PWM  duty  cycle  inputs. 


Figure  4.  Signal  Flow  (Transmitter  to  Servo  Output) 


The  servo  rotations  are  then  transmitted  through  linkages  that  control  the 
helicopter  control  surfaces.  More  specifically,  the  main  rotor  linkages  come  from  the 
cyclic  (lateral  and  longitudinal)  servos  and  from  the  collective  servo.  Figure  5  illustrates 
the  signal  flow  Ifom  the  servo  input  to  the  commanded  flap  angle  as  well  as  the  signal 
flow  from  the  commanded  flap  angle  to  the  actual  flap  angle.  The  linear  flapping 
dynamics  are  explained  in  3.B.7. 


1-45'  145] 


Linkage 


Pi 

[-20'  +2(7] 


Linear  Flapping 
Dynamics 


Figure  5.  Signal  Flow  (Servo  to  Commanded  Flap  Angle) 
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The  linkage  that  transforms  the  servo  mechanical  rotation  to  a  control  surface 
command  is  linear  and  proportional  to  the  servo  mechanical  rotation. 

Figure  6  illustrates  the  model,  which  simulates  the  channel,  from  control  input 
m  =  ]  to  PWM  duty  cycle  with  limits  =  [600  1520  2400] 


Figure  6.  Channel  Model. 


The  first  four  signals  in  Figure  6  represent  the  command  inputs  u  ;  the  last  several 
signals  are  other  signals  that  are  being  forwarded  to  the  model.  Saturation  limits  are 
placed  according  to  the  control  input  limits.  The  .5  constant  on  the  third  signal  ensures 
that  when  the  control  effort  is  zero,  the  UAV  will  maintain  a  50%  collective  input  for 
hovering  flight.  This  is  because  the  main  rotor  collective  is  designed  for  hovering  at  50% 
command  input.  Figure  7  illustrates  the  model  of  the  helicopter  servos  with  PWM  duty 
cycle  as  inputs  and  mechanical  servo  rotations  as  outputs. 
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Figure  7.  Servo  Model  (Control  Input  to  PWM) 


Each  servo  block  also  has  second  order  filters  that  model  the  time  responses  of  the 
mechanical  output.  The  time  response  models  for  the  servos  were  obtained  from 
[1,  p.  52],  The  time  response  model  for  the  servo  mechanical  output  is  illustrated  in 
Figure  8. 


( 1 : 

pwm 


Figure  8.  Servo  Input  to  Output  Filter 

Figure  9  illustrates  the  PWM,  servo  mechanical  rotation,  and  control  surface 
command  outputs  from  several  longitudinal  step  commands  of  -1,  0,  and  1.  These  inputs 
cover  the  full  range  of  stick  input  motion  for  longitudinal  control. 
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Figure  9.  Channel  Model  Responses 

Initially,  when  the  input  is  zero,  the  duty  cycle  output  is  in  the  middle  of  the  duty 
cycle  range,  1520  /us ,  which  gives  the  zero  position  on  the  servo  and  no  control  surface 
command.  When  the  longitudinal  input  is  maxed  out  at  full  forward  (+1),  the  duty  cycle 
maxes  out  at  2400  jus  and  this  generates  the  maximum  servo  rotation  of  -t-45°  and  a 
proportional  cyclic  control  surface  command  of  20°.  It  is  evident  that  the  full  back 
command  (-1)  generates  the  -20°  collective  control  surface  command.  The  second  order 
filter  models  the  second  order  behavior  of  the  servo  mechanism  and  gives  a  settling  time 
of  .055  seconds  with  7.6%  overshoot. 
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The  duty  cycle  limits,  servo  mechanical  outputs,  and  other  parameters  that  are 
accessed  by  the  table  blocksets  in  the  models  (illustrated  in  Figures  5  through  7)  are 
defined  in  “Config.m”  outlined  in  Appendix  B.  This  file  also  defines  other  constant 
parameters  like  the  UAV  mass,  dimensions,  moments  of  inertia,  and  so  on. 

B,  MAIN  ROTOR  DESIGN 

The  main  rotor  thrust  is  directly  related  to  the  pilot  commanded  collective  and 
indirectly  related  to  the  vehicle  velocity  states  (w,  v,  w)  and  the  rotor  inflow  velocity 
( V. ).  This  relation  is  given  by  Padfield  [4]  and  Leishman  [5]  as; 


T  = 


0^ 


2  ^ 


Az+^, 


(3.1) 


Determination  of  the  thrust  is  complicated  by  the  fact  that  the  rotor  inflow,  ,  is  directly 
influenced  by  the  thrust.  Padfield  explains  that  as  the  rotor  thrust  increases,  the  rotor 
inflow  also  increases  [4,  p.  122].  The  opposite  is  also  true.  That  is,  as  the  rotor  thrust 
decreases,  the  main  rotor  inflow  velocity  decreases.  Knowing  this  helps  the  designer  to 
understand  how  the  thrust  and  inflow  interact  and  the  ultimate  behavior  of  the  thrust  to  a 
pilot  commanded  step  in  collective,  which  is  denoted  by  0^ . 

If,  for  example,  the  pilot  gives  a  step  collective  command,  increasing  ,  the  thrust 
will  begin  to  increase.  This  occurs  because  the  rotor  inflow  term,  X. ,  has  not  had  enough 

time  to  increase  to  the  critical  value  that  will  cause  the  thrust  to  decrease  to  a  new 
equilibrium.  So,  during  this  phase  of  the  thrust  response,  the  thrust  is  increasing  since  the 
change  in  0^  is  larger  than  any  increase  in  /I. ;  note  that  the  X.  term  is  being  subtracted  in 

Equation  (3.1).  Later,  it  will  be  shown  how  the  rotor  inflow  increases  with  increasing 
thrust  but,  for  now,  the  thrust  response  to  an  increasing  inflow  will  be  analyzed.  As  the 
rotor  inflow  velocity  increases,  the  X.  term  gets  larger.  The  same  is  true  for  the  term 

when  the  vehicle  is  climbing.  This  causes  the  thrust  to  begin  to  decrease  once  X.  reaches 
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a  critical  value.  The  thrust  will  continue  to  deerease  until  it  reaches  a  new  equilibrium 
point.  To  better  understand  this  behavior,  the  process  for  determining  the  rotor  inflow 
will  now  be  discussed. 

1.  Rotor  Inflow  from  Momentum  Theory 

Momentum  analysis  in  forward  flight  provides  for  the  solution  of  the  main  rotor 
inflow  in  forward  flight  as  well  as  in  elimb.  Due  to  the  eomplexity  of  the  main  rotor 
inflow  in  descent,  a  different  technique  will  be  used  for  this  flight  configuration. 

The  relationship  between  rotor  thrust  and  the  rotor  inflow  is  dietated  by  the 
conservation  of  mass  and  the  conservation  of  momentum  and  energy  [4,  pp.  115-116]. 
The  theory  that  allows  for  the  derivation  of  these  relationships  is  known  as  momentum 
theory.  Conservation  of  mass  dictates  that  the  mass  flow  rate  of  a  volume  of  air  passing 
through  the  main  rotor  disk  area  is  proportional  to  the  volume  flow  rate  (m^/sec)  and  the 
air  density.  This  relation  is  given  by  Anderson  as  [6,  p.  93]; 

m  =  p(V„A,)  (3.2) 

where  is  the  normal  component  of  the  velocity  of  the  mass  of  air  through  the  rotor 

disk.  The  thrust  generated  from  the  aceeleration  of  the  air  through  the  main  rotor  is 
simply 

T  =  ma  (3.3) 

Padfield  [4]  relates  the  rate  of  ehange  of  momentum  between  the  undisturbed  upstream 
eonditions  and  the  far  wake  to  the  rotor  loading  and  is  given  by  [4,  p.  116]; 

T  =  mv.^  (3.4) 

where  is  the  induced  flow  far  downstream  of  the  rotor.  Padfield  proeeeds  to  relate  the 

ehange  in  kinetic  energy  of  the  flow  to  the  work  done  by  the  rotor,  leading  to  the 
eonclusion  that  the  velocity  of  the  air  far  downstream  of  the  rotor  is  twice  the  induced 
velocity  of  the  air  immediately  upstream  of  the  rotor  blades.  This  gives  the  expression 
[4,p.  116]; 

fi.=2v,  (3.5) 
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Substituting  Equations  (3.2)  and  (3.5)  into  (3.4)  gives  [4,  p.  116]; 


T  =  (pA,FJx(2v,)  (3.6) 

The  veloeity  eomponent  normal  to  the  rotor  depends  on  the  helieopter  flight 
configuration.  In  the  hover  flight  configuration,  the  upstream  velocity  of  the  air  normal  to 
the  rotor  is  the  same  as  the  inflow  velocity,  therefore,  =  v.  in  a  hovering  flight 
configuration.  In  a  climbing  flight  configuration,  the  upstream  velocity  of  the  air  normal 
to  the  rotor  is  +  v- .  In  a  descent,  the  upstream  velocity  of  the  air  normal  to  the  rotor  is 
-V.  [4,  p.  116].  At  low  forward  airspeeds,  the  normal  component  is  given  by  Prouty  as 

[4,  p.  123].  Higher  airspeed  forward  flight  is  the  more  general  case  and  simpler 

to  use  because  the  normal  component  velocity  reduces  to  +  v]  at  the  lower  airspeeds. 

The  velocity  component  normal  to  the  disk  rotor  for  all  forward  airspeeds  is  given  by 
Leishman  as  [5,  p.  63]; 

{Hovering  and  Forward  Airspeed) 

r-2 - ^ ^  (3.7) 

v„=^jv^  +2Vv.sm  a +vf 

Where  a  is  the  angle  between  the  horizon  and  the  main  rotor  blade  tip  path  plane. 
Equation  (3.7)  is  useful  because  it  reduces  to  the  normal  velocity  component  for  hovering 
flight  as  well  as  forward  flight  when  the  helicopter  is  in  either  of  these  configurations. 
Rearranging  (3.6)  for  v.  gives; 


T 


(3.8) 


When  hovering,  the  velocity  of  the  air  normal  to  the  rotor  ( T„ )  is  equal  to  the  rotor  inflow 
(v.).  Substituting  y.  for  in  (3.8)  gives  T  =  ihover'^P^d  ■  Substituting  this  into  (3.8) 
gives; 
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T  _  ^ihover  )  _  ^ihover 

2pA,V„  (2pA,)V„  ■  F„ 


(3.9) 


Substituting  (3.7)  into  (3.9)  gives  the  main  rotor  inflow  for  any  airspeed  and  orientation 
as  [5,  p.64]; 


V,.  = 


2 


V 


ihover 


COS 


(«))  +(^"ooSin(a)  +  v.) 


(3.10) 


The  velocity  parallel  to  the  main  rotor  plane  is  cos(«) .  Normalized  by  the  main  rotor 
tip  speed  (Qi?)  gives  the  normalized  translational  velocity  // =  cos(a)/(Q7?) .  The 

main  rotor  inflow  ratio  is  [5,  p.  65]; 


^  V, 

OR  OR 


(3.11) 


Leishman  then  writes  the  first  term  in  (3.11)  in  terms  of  the  normalized  normal  airspeed 
p  by  relating  sin(«)  =  cos(a)tan(«)  = //tan(«) .  Therefore,  Equation  (3.11)  is 
now 

/I  =  //tan(a)  +  -^  (3.12) 

FIR 


Substituting  (3.10)  into  (3.11)  gives 


1 

;i  =  //tan(a)  +  -pJi^^  (3.13) 

The  inflow  ratio  at  hover  is  given  by  [4]; 


The  solution  to  the  inflow  ratio  is  then  [4,  p.  65] 
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(3.14) 


A,  =  //tan(«)  + 


a 


Leishman  explains  that  a  numerical  solution  to  Equation  (3.14)  can  be 
implemented  using  the  Newton-Raphson  technique.  The  advantage  to  using  this  method 
is  a  more  rapid  convergence  as  well  as  the  fact  that  it  can  be  employed  for  practical 
calculations  involving  rotors  in  all  ranges  of  flight  configurations,  from  hover  to  climb,  as 
well  as  forward  flight  and  descents.  Equation  (3.14)  does  provide  convergence  in 
descents;  however,  the  solution  may  not  be  physically  realizable.  In  descent  cases,  the 
rotor  inflow  velocity  will  be  approximated  as  linear  and  will  be  discussed  in  Chapter 
III.B.4. 


2,  Newton-Raphson  Iteration  Technique 


The  Newton-Raphson  iteration  scheme  is  implemented  in  the  estimation  of  the 
main  rotor  and  tail  rotor  inflow.  The  inflow  is  iteratively  estimated  from  the  previous 
value  until  the  error  between  the  current  and  previous  value  is  less  than  some  acceptable 
tolerance.  The  inflow  is  determined  iteratively  using  Equations  (3.15)  through  (3.18) 

[5,p.  67]; 


/W 

k'wi 


(3.15) 


f{X)  =  X-/uiMi[a) - ,  (3.16) 

/(i)  =  l  +  ^(y+i^pi  (3.17) 


The  iteration  scheme  continues  until  the  error  estimator  is  less  than  .005%.  The 
error  estimator  is  given  as  [5,  p.  66]: 


s  = 


^n+\  \ 


A 


«+l 


(3.18) 
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For  any  given  thrust,  Equation  (3.15)  will  converge  to  the  correct  value  of  inflow  by 
using  an  initial  inflow  estimate  of  \  .  The  hover  inflow  is  determined  at 

hovering  conditions,  where  the  main  rotor  thrust  is  known  to  be  equal  to  the  weight  of  the 
helicopter.  Since  =  v.  at  hover.  Equation  (3.8)  can  be  rewritten  as: 


V,.  =  ■ 


v!=- 


2pA, 


*^ihover 


'2pA, 


(3.19) 


where  the  weight  of  the  helicopter  is  W  =  8.2  x  9.806  =  80.41  N  .  The  main  rotor  disk  area 
is  A^  =  =  7rx(^.775y  =1.887  ,  which  gives  a  main  rotor  hovering  inflow  of 

V.  =  4.171  m/s  .  Normalized  by  the  rotor  tip  speed  (Qi?)  gives  =32.23x10  ^ 

The  following  example  shows  the  inflow  convergence  for  both  cases  where  the 
thrust  is  increasing  and  decreasing.  Initially,  the  main  rotor  thrust  increases  with  an 
increase  in  collective  lever  input.  Eater,  it  will  be  shown  how  this  occurs,  but  for  now  it 
will  be  assumed  that  the  main  rotor  thrust  for  two  consecutive  samples  are  T=[80.5  8E3]. 
There  is  no  translational  velocity  component,  so  p  =  0  and  the  tip  path  plane  angle, 
a  =  0  .  The  coefficients  of  thrust  for  these  two  samples  are  obtained  by  normalizing  the 
thrust  values  of  each  sample  by  p{QRf  and  gives 

Q  =  [^2.079x10^^  2.0997x10^^^.  Running  these  thrust  values  through  “getvi.m” 

outlined  in  Appendix  B,  gives  the  two  convergence  trajectories  illustrated  in  Eigures  9 
and  10.  The  m  files  used  for  generating  convergent  inflow  velocities  for  the  main  rotor 
and  tail  rotor  are  listed  in  Appendix  B. 
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Figure  10.  Convergent  Inflow  Using  Newton-Raphson  Technique  (T=80.5  N) 
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Figure  1 1 .  Convergent  Inflow  Using  Newton-Raphson  Technique  (T=8 1 .3  N) 


Figure  10  illustrates  the  value  of  the  estimated  rotor  inflow  for  up  to  6  iterations. 
The  inflow  estimate  starts  at  =4.171  m/s,  before  the  first  iteration.  The  Newton- 

Raphson  iteration  loop  in  “getvi.m”  increases  the  value  of  v.,  using  Equations  (3.15) 
through  (3.17),  until  s  <  .005% .  As  soon  as  the  error  becomes  less  than  .005%,  the  value 
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of  the  inflow  for  that  iteration  routine  is  extraeted  and  used  to  eompute  the  new  value  of 
thrust.  From  this  example,  it  is  elear  that  as  the  thrust  inereases  the  main  rotor  inflow 
veloeity  also  increases.  The  main  rotor  inflow  increases  from  the  initial  estimate  of  4.171 
m/s  to  4.175  m/s,  which  is  the  new  rotor  inflow  at  T=80.5  N.  As  the  thrust  increases  to 
81.3  N  on  the  second  sample,  the  rotor  inflow  converges  to  4.215  m/s,  as  illustrated  in 
Figure  11.  Figure  12  illustrates  the  main  rotor  inflow  convergence  when  the  rotor  thrust 
decreases. 
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Figure  12.  Convergent  Inflow  Using  Newton-Raphson  Technique  (T=79.0  N) 


In  this  case,  as  the  thrust  decreases  to  79.0  N,  the  main  rotor  inflow  velocity 
decreases  to  4.0971  m/s.  This  iteration  scheme  is  implemented  in  the  Simulink  model  for 
both  the  main  rotor  and  tail  rotor  at  a  rate  of  100  Hz. 


3,  Main  Rotor  Thrust  and  Inflow  Response 

The  helicopter  simulation  starts  in  a  hovering  configuration  since  the  inflow  and 
thrust  values  are  known  in  this  regime.  As  the  main  rotor  collective  is  changed,  the  rotor 
inflow  is  determined  as  described  in  section  III.B.2.  The  convergent  inflow  is  then  used 
in  Equation  (3.1)  to  calculate  the  main  rotor  thrust. 
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As  the  pilot/operator  gives  a  collective  input  command  {0^),  the  new  value  of 


thrust  is  determined  from  Equation  (3.1)  by  using  the  previous  known  value  of  the  inflow 
(1  sample  delay).  Then,  the  new  value  of  thrust  is  input  into  the  Newton-Raphson 
technique  to  generate  the  new  value  of  the  inflow  and  the  process  starts  over  again;  the 
loop  is  closed  in  this  fashion. 

As  explained  earlier,  Equation  (3.1)  reveals  that  as  the  inflow  increases  in 
response  to  an  increase  in  thrust,  the  thrust  will  begin  to  decrease  because  of  the 

f  u  +;i'\ 

— - -  term  that  is  being  subtracted.  A  simulation  was  run  in  order  to  illustrate  this 

2  j 

behavior.  The  input  is  a  step  collective  from  5.495°  to  5.657°.  The  response  of  the  thrust 
and  main  rotor  inflow  velocity  is  illustrated  in  Eigure  13. 
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Rotor  Inflow  Response 
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Eigure  13.  Main  Rotor  Thrust  and  Inflow  Response  to  Step  Collective 


The  thrust  response  illustrated  in  Eigure  13  shows  that  the  thrust  returns  to  the 
equilibrium  thrust.  This  response  is  in  agreement  with  [4,  p.  305],  which  illustrates  the 
thrust  response  to  a  step  collective  for  a  full-scale  helicopter.  In  this  case,  the  equilibrium 
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value  of  the  thrust  is  the  weight  of  the  vehicle.  This  behavior  guarantees  the  first  order 
heave  response  that  is  characteristic  of  both  full-scale  and  miniature  helicopters.  That  is, 
for  any  given  step  collective  input,  the  velocity  rate  in  the  z  direction  will  display  a 
second-order,  over  damped  response.  Figure  14  illustrates  the  heave  rate  response  for  a 
collective  step  input  from  5.495°  to  6.5° .  This  heave  response  is  also  in  agreement  with 
the  vertical  axis  response  characteristic  detailed  in  [4,  p.  401]. 
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Figure  14.  Fleave  Response  to  Collective  Command  Simulation 

Figure  14  illustrates  that  since  the  thrust  returns  to  Thrust=80.41  N,  after  the  peak, 
there  is  no  net  acceleration,  and  this  is  why  the  velocity  settles  to  a  constant  value. 
Equation  (3.1)  in  the  helicopter  model  is  implemented  in  “UpdateMR.m,”  outlined  in 
Appendix  B. 
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4,  Main  Rotor  Thrust  and  Inflow  Response  in  Pure  Descent 

There  exist  several  operating  conditions  within  the  descending  flight 
configuration.  These  are  low  rate  descent,  Vortex  Ring  State  (VRS),  and  the  Windmill 
Break  State  (WBS)  (also  known  as  autorotation).  In  low  rate  of  descent,  there  is  a  well- 
defined  slipstream  in  the  downward  direction  and,  in  this  case,  momentum  theory  can  be 
applied.  In  this  flight  regime,  the  induced  velocity  is  slightly  increasing  with  increasing 
descent  velocity.  This  region  of  allowable  low  rate  descent,  in  the  absence  of  forward 
airspeed,  is  small.  In  the  WBS,  the  slipstream  is  in  the  upward  direction  and  is  a 
condition  where  the  energy  imparted  onto  the  rotors  by  the  slipstream  is  sufficient  to 
drive  the  rotor.  The  slipstream  is  well  defined  and,  momentum  theory  can  also  be  applied 
in  this  state. 

Once  the  descent  velocity  breaks  the  limit  afforded  by  momentum  theory,  the 
vehicle  enters  what  is  known  as  the  Vortex  Ring  State  (VRS).  Padfield  explains  that 
[4,p.  102]; 

Flight  test  experience  on  a  small  tandem  rotor  helicopter... reported  that 
the  characteristic  vibration  [of  VRS]  began  at  a  rate  of  descent  equal  to 
about  23%  of  the  hovering  induced  velocity  and  persisted  until  the  rate  of 
descent  exceeded  125%  of  the  hover-induced  velocity. 

This  boundary  cannot  be  known  a-priori  and  can  only  be  established  during  flight- testing. 
This  being  said,  the  model  will  not  be  able  to  accurately  simulate  vehicle  dynamics  once 
the  VRS  is  entered  into,  and  will  therefore  be  avoided  during  simulations.  Padfield’s 
referenced  flight  test  is  used  to  set  an  initial  VRS  boundary  for  the  X-Cell  model.  A  more 
conservative  descent  limit  of  20%  hovering  induced  velocity  gives  a  VRS  upper 
boundary  of  =  .834  m/s .  The  lower  VRS  boundary  is  =  5.2  m/s  .  Model  simulation 

is  not  concerned  with  aerobatic  maneuvers;  therefore,  a  pure  descent  command  will  not 
exceed  .834  m/s.  It  should  be  mentioned  at  this  point  that  VRS  can  be  avoided  entirely, 
within  the  VRS  descent  boundary,  by  providing  supplemental  forward  airspeed 
commands. 

In  order  to  generate  a  smooth  simulation,  Padfield’s  linear  approximation  is  used 
in  the  VRS  region.  These  are  given  as  [4,  p.  118],  [5,  p.  57]; 
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Equation  (3.20)  describes  the  main  rotor  inflow  when  the  vehicle  is  within  the  upper 
descent  boundary,  between  hover  and  -.834  m/s.  Equation  (3.21)  describes  the  main  rotor 
inflow  when  the  vehicle  is  descending  at  a  rate  between  -.834  m/s  and  -8.34  m/s. 
Equation  (3.22)  describes  the  main  rotor  inflow  when  the  vehicle  is  descending  faster 
than  -8.34  m/s. 


5,  Main  Rotor  Inflow  with  Increasing  Airspeed 

Eigure  15  illustrates  that  for  small  tip  path  plane  angles  (a  =  O) ,  where  the  drag  is 

minimal,  the  convergent  inflow  velocity  decreases  with  increased  airspeed.  As  the  drag 
becomes  more  and  more  dominant,  as  is  the  case  with  «  =  6° ,  the  inflow  begins  to 
increase  again  at  some  critical  airspeed. 


Main  Rotor  Inflow  Velocity  vs.  Airspeed:  Forward  Flight 


Eigure  15.  Convergent  Inflow  Velocities  in  Eorward  Elight  (After  [5]) 
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In  most  applications,  the  tip  path  angle  of  the  main  rotor  will  not  be  6°  when  the 
vehicle  is  hovering  (zero  airspeed).  In  order  to  compensate  for  drag,  the  tip  path  angle 
will  only  increase  steadily  as  the  vehicle  airspeed  increases.  Figure  16  is  a  better 
depiction  of  the  convergent  rotor  inflow  velocity  behavior  for  the  model,  over  the  range 
of  operational  airspeeds. 
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X-Cell  60  Main  Rotor  Inflow  Velocity  vs.  Airspeed 
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Figure  16.  Convergent  Inflow  Velocities  of  X-Cel  60  (After  [5]) 


Figure  16  takes  into  account  that  the  drag  and  tip  path  angle  increases  as  the 
airspeed  increases.  Figure  16  is  labeled  in  several  locations  to  show  the  drag  and  tip  path 
plane  angle  at  that  specific  operating  condition.  For  example,  at  a  forward  airspeed  of  1 
m/s,  the  drag  is  approximately  .1  Newtons  with  a  tip  path  plane  angle  of  .07°  and  a 
convergent  inflow  velocity  of  4.11  m/s.  At  a  forward  airspeed  of  5  m/s,  the  drag  is  2.4 
Newtons,  the  tip  path  plane  angle  is  1.7°  with  an  inflow  velocity  of  3.11  m/s.  It  can  also 
be  seen  that  at  approximately  8  m/s,  the  inflow  begins  to  increase  again  and  gets  larger 
than  the  hover  inflow  velocity  at  around  13.5  m/s  (30.2  mi/hr). 

Figure  16  also  represents  the  range  of  main  rotor  inflow  velocities  for  the  vehicle 
in  trimmed  flight.  That  is  to  say,  that  the  vehicle  will  generate  an  inflow  of  3.1 1  m/s,  for 
example,  at  a  steady  state  trimmed  flight  of  5  m/s  forward  airspeed.  As  the  vehicle 
accelerates  to  different  operating  points,  there  will  be  inflow  transients  in  between  the 
operating  points. 
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6.  Main  Rotor  Torque 

At  hover,  the  coefficient  of  torque  is  a  combination  of  the  induced  torque 
coefficient  and  the  profile  torque  coefficient.  The  induced  torque  coefficient  is  due  to  the 
[5]; 

. .  .combined  effect  of  the  rearward  tilt  of  the  incremental  lift  vectors. 


This  is  a  type  of  aerodynamic  drag  on  the  blades,  which  is  proportional  to  the 
thrust  and  is  expressed  as  [3,  p.  22] 


(3.23) 


The  profile  torque  can  be  thought  of  as  coming  from  friction  on  the  blade  surface  and  is 
expressed  as 
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Prouty  presents  where  the  drag  coefficient  is  treated  as  a  constant, 

c"”'  =  .024 

d 

Combining  both  expressions  gives  the  total  torque  coefficient  as  [3,  p.  22]; 
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In  forward  flight,  the  main  rotor  torque  coefficient  is  [3,  p.  133]; 
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The  yaw  moment  generated  by  this  torque  is  given  as  [3,  p.  22]; 
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(3.27) 


A'„=p(aR)’^Cf 

The  following  is  a  treatment  of  the  main  rotor  yaw  moments  that  are  generated  at 
various  operating  conditions; 

Yaw  Moment  at  3  m/s  Forward  Flight: 

r  =  81.1725  jV 


Yaw  Moment  at  Hover: 


r  =  80.4145  jV 


p{aR)  A, 


=.0020769 


=  2.0986x10^ 
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=.0020964 


=  2.1081x10 
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1  jyZ 

=  p{nRf  -^—Cq  =  6.326  Nm 


The  conclusion  from  this  exercise  is  that  the  forward  airspeed  is  not  going  to  have 
significant  change  on  the  rotor  torque.  The  general  trend  is  that  the  yaw  moment  will 
increase  as  the  main  rotor  thrust  increases. 


7.  Main  Rotor  Flapping  Dynamics 

The  main  rotor  flapping  dynamics  displays  a  second  order,  over  damped  response. 
Equation  (3.28)  expresses  that  the  main  rotor  flapping  is  determined  by  the  longitudinal 
and  lateral  flap  angle  commands,  and  ,  respectively.  These  commanded  flap 
angles  are  proportional  to  the  commanded  inputs  and  u,^,.  The  direct  relation 
between  u  and  ^  is  discussed  in  Chapter  III.A.  Additionally,  roll  and  pitch  motion  from 
the  UAV  affects  the  main  rotor  flapping  directly.  Mettler  gives  the  linear  flapping 
dynamics  as  [2,  p.  81]; 
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where  is  the  main  rotor  time  constant  to  be  determined  in  system  identification.  The 

value  of  =  .l  seconds  was  obtained  from  [2]  at  hovering  conditions.  The  eigenvalues 

of  the  state  transition  matrix  A  in  (3.28)  gives  two  stable  eigenvalues/poles  at  -10.  Taking 
the  A  and  B  matrices  and  defining  the  output  C  matrix  as  C=[l  0]  with  input  coupling 
D=[0  0  0  0],  the  SISO  transfer  function  for  this  linear  system  is  determined  using  the 
following  Matlab  commands; 

»  [num,den]=ss2tf(a,b,c,d,3) 

»sys=tf(num,  den) 

The  3  in  the  argument  of  the  first  function  defines  the  input  as  ,  which  is  the 
third  input  for  a  four  input  system.  The  transfer  function  representing  Equation  (3.28) 
from  step  command  (5'" to  flap  angle  is  then  given  as; 

Plan  _  lO.y  +  100  P29) 

/ +205  +  100  ■  ^ 

An  approximation  to  the  settling  time  for  this  second  order  system  is 

7;=^  (3.30) 

where  the  denominator  of  Equation  (3.29)  gives 

+2<^Cf)„+col  +20s  +  m  (3.31) 

Solving  for  the  natural  frequency  and  damping  ratio  gives  settling  time  estimate  of 
7)  =  .4  seconds,  which  is  in  agreement  with  the  responses  illustrated  in  Eigures  16  and  17. 

Eigures  16  and  17  illustrate  the  main  rotor  flapping  responses  to  longitudinal  and 
lateral  input  commands  as  well  as  the  vehicle  velocity  responses  in  the  lateral  and 
longitudinal  directions.  The  lateral  step  response  of  the  velocity  in  Eigure  17  makes  sense 
for  the  first  2.5  seconds.  The  drop  in  velocity  after  2.5  seconds  is  due  to  the  tendency  of 
the  tail  boom  to  rotate  away  from  the  direction  of  lateral  turn,  in  the  absence  of  a  tail 
rotor  correction  command. 
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Step  Command  vs.  Flap  Angle  Response  (deg) 
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Figure  17.  Lateral  Flap  Angle  Response 
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Figure  18.  Longitudinal  Flap  Angle  Response 
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The  determination  of  the  velocity  responses  illustrated  in  Figures  16  and  17 
originate  from  the  thrust  vector  being  translated  into  the  x  and  y  axes  as  the  main  rotor 
flaps.  This  will  be  discussed  in  the  next  section. 


8. 


Main  Rotor  Forces  and  Moments 


In  order  to  better  appreciate  and  understand  the  forces  and  moments  acting  on  a 
helicopter,  a  free  body  diagram  is  illustrated  in  Figure  19. 


Figure  19.  Main  Rotor  Forces  and  Moments  (After  [2,  p.  75]) 


Figure  19  illustrates  the  forces  and  moments  acting  on  the  helicopter  as  the  main 
rotor  tilts  from  a  flap  angle  command.  As  the  operator  inputs  a  forward  cyclic  command, 
the  helicopter  main  rotor  will  provide  a  forward  tilt,  which  generates  the  longitudinal  flap 
angle  illustrated  as  (-a)  in  Figure  19.  As  the  main  rotor  tilts  forward,  the  main 

rotor  thrust  vector  also  tilts,  thereby  generating  a  force  in  the  +x  direction  (7)^^) . 
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The  thrust  (J)^^^)then  causes  an  acceleration  in  the  +x  direction.  Additionally,  the  x 

component  of  the  thrust  vector  acts  on  the  moment  arm  between  the  main  rotor  hub  and 
the  vehicle  center  of  gravity,  which  in  turn  generates  the  aerodynamic  pitching  moment 

As  the  main  rotor  tilts,  there  is  a  restoring  moment  due  to  the  centrifugal  forces 
generated  by  the  rotating  blades.  For  now,  the  restoring  moment  that  is  generated  from 
the  main  rotor  centrifugal  force  is  modeled  as  a  spring  moment  (M j^);  this  spring 

moment  is  modeled  as  being  proportional  to  the  flap  angle  .  An  additional  moment, 

not  illustrated  in  Figure  19,  is  a  damping  moment  which  is  proportional  to  any  pitching  or 
rolling  rate,  q  ox  p  .  This  damping  moment  is  physically  provided  by  the  stabilizer  bar, 
which  behaves  as  a  main  rotor  cyclic  control  augmentation  device.  The  main  rotor 
stabilizer  bar  does  not  generate  any  lift;  instead,  it  is  a  type  of  gyroscope  in  the  sense  that 
it  maintains  a  horizontal  orientation  given  any  main  rotor  orientation.  As  the  main  rotor 
flaps  in  any  direction,  the  stabilizer  bar  maintains  a  horizontal  attitude.  Linkages  between 
the  stabilizer  bar  and  main  rotor  control  surfaces  cause  the  augmented  control  inputs  from 
the  stabilizer  bar  to  correct  for  any  roll  or  pitch  rates  in  order  to  drive  the  rates  to  zero. 
The  pitch  damping  moment,  for  example,  is  proportional  to  the  pitch  rate  q  and  is 
expressed  as  [2]; 

(3.32) 

where  and  are  constants  to  be  determined  by  system  identification.  The  term 

describes  the  rate  at  which  the  vehicle  pitches  given  a  perturbation  in  the  flap  angle  . 

See  Appendix  A  for  identification  model  parameters  used  from  the  MIT  Instrumented  X- 
Cell  60.  The  X-Cell  60  model  parameters  at  hover  are  used  to  model  the  unknown  T-Rex 
Align  parameters  until  system  identification  is  completed. 

Conducting  a  force  and  moment  analysis  for  all  6-DOF  gives  Equations  (3.33) 
and  (3.34),  for  the  summation  of  all  forces  and  moments. 
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(3.33) 


^  MR  ~  '^MR  sin(Ao«)cos(^,,,) 

^MR  ~  '^MR  sin(y5to)cos(^,„„) 

^  MR  ~  "^MR  cos(A„„)cos(^,,,) 

^MR  ~  ^MR^MR  +  i^p)Plat  ~  J  P 
^ MR  ~^MR^MR  '^i^p)Plon  ~^a'^fP  (3.34) 

Nmr-p{^)"^Cq 

It  should  be  noted  that  the  Simulink  model  implementation  of  Equation  (3.34) 
takes  into  eonsideration  the  sign  of  the  flap  angle  before  applying  the  equation.  That  is,  if 
the  longitudinal  flap  angle  is  positive,  the  restorative  moment  {Kp)xp,^^  must  be 

negative;  if  the  flap  angle  is  negative,  the  restorative  moment  {t^p)'^Pion  must  be 

positive.  This  ensures  that  the  aerodynamie  pitehing  moment  and  the  restorative  moments 
are  opposite  in  sign  and  will  eaneel  eaeh  other.  Note  that  the  yaw  axis  main  rotor  moment 
(^Affl)  was  eovered  in  Chapter  III.B.6. 

9,  Pitch/Roll  Rate  and  Velocity  Response  to  Main  Rotor  Forces  and 
Moments 

Figure  20  illustrates  the  piteh  and  veloeity  response  generated  by  applying 
equations  (3.33)  and  (3.34)  to  a  lateral  right-stiek  eyelie  input.  These  responses  are  in 
elose  agreement  with  the  full-seale  helieopter  response  in  [4,  p.  349]. 
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Percent  Lateral  Cyclic  Step  Input 


Figure  20.  Lateral  Cyelie  Step  Response 


Figure  20  illustrates  that  a  step  input  in  lateral  eyelie  generates  an  inereasing  roll 
rate.  The  differenee  between  the  model  roll  rate  response  and  the  full  seale  helieopter 
response  of  [4]  is  that  the  model  response  is  more  stable  in  that  the  roll  rate  damps  baek 
to  zero  rate  while  the  eyelie  input  is  still  aetive.  This  is  highly  desirable  and  eharaeteristie 
of  the  miniature  helieopter  model.  If  the  roll  rate  is  too  damped,  the  roll  derivative  ean 

be  tuned,  further,  this  will  be  determined  in  future  system  identifieation  of  the  T-Rex 
Align.  Figure  20  also  illustrates  that  the  roll  angle  only  attains  a  maximum  of  .  2°  roll 
with  a  damping  baek  to  zero  roll  angle.  This  is  also  highly  desirable  beeause  it  allows  the 
vehiele  to  eontinue  to  inerease  its  veloeity  in  the  y  direetion  without  being  hindered  by 
inereasingly  large  roll  angles  generated  by  the  rolling  moment.  This  illustrates  that  the 
restorative  and  damping  moments  do  a  good  job  of  eaneeling  aerodynamie  moments 
generated  by  the  flapping  dynamies.  Figure  20  also  shows  that  the  veloeity  in  the  y 
direetion  inereases  as  the  eyelie  input  is  maintained  at  the  20%  position.  The  longitudinal 
responses  are  similar  to  those  illustrated  in  Figure  20. 
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C.  TAIL  ROTOR  DESIGN 

The  tail  rotor  behaves  mueh  the  same  way  the  main  rotor  does.  Therefore,  the  tail 
rotor  thrust  and  inflow  are  first  determined  in  the  hover  conditions.  A  tail  rotor  collective 
command  will  then  cause  the  tail  rotor  thrust  to  either  increase  or  decrease,  depending  on 
the  direction  of  the  command.  The  rotor  inflow  will  either  increase  with  increased  tail 
rotor  trust  or  decrease  with  decreased  tail  rotor  thrust.  As  the  tail  rotor  inflow  velocity 
increases  with  increased  thrust,  for  example,  the  tail  rotor  thrust  will  first  peak  then  it  will 
begin  to  decrease  just  as  was  the  case  with  the  main  rotor.  The  same  iterative  procedure 
will  be  applied  to  determine  the  convergent  inflow  for  the  tail  rotor.  There  are  several 
other  differences,  however,  that  need  to  be  discussed. 

1,  Tail  Rotor  Hub  Airspeeds  and  Main  Rotor  Wake  Considerations 

Modeling  of  the  tail  rotor  requires  an  understanding  of  how  the  helicopter’s 
motion  affects  the  local  tail  rotor  airspeed  components.  In  pure  translational  motion,  the 
tail  rotor  hub  airspeed  is  straightforward  to  calculate.  If  the  vehicle  is  flying  in  a  pure 
sideways  motion  (positive  y  -axis),  the  velocity  of  the  tail  will  be  the  same  as  the  vehicle 
sideways  airspeed.  If  the  vehicle  is  experiencing  roll,  pitch,  or  yaw  rates,  this  will  affect 
the  tail  rotor  velocity  directly.  The  tail  rotor  hub  airspeed  in  the  y  direction  is  given  by 
Gavrilets  as  [1,  p.  51]; 


v,r^v^-hrr-K\p\  (3.35) 

where  is  the  resultant  airspeed  of  the  vehicle  in  the  y  direction,  compensated  for  wind 
disturbance.  Equation  (3.35)  expresses  the  tail  rotor  hub  velocity  as  a  combination  of  the 
translational  velocity  of  the  tail  (v^),  angular  velocity  of  the  tail  due  to  tail  rotor  yaw 
( l,^r ),  and  angular  velocity  of  the  tail  due  to  roll  ( ).  If  the  vehicle  travels  in  the 
positive  y  direction,  a  positive  yaw  rate  (clockwise)  will  reduce  the  net  tail  rotor  airspeed 
in  the  y  direction.  If  the  yaw  rate  is  negative  (counterclockwise,  with  -r  ),  the  tail  rotor 
airspeed  becomes  meaning  that  the  tail  rotor  hub  airspeed  is  faster  than  the 

vehicle  center  of  gravity.  The  tail  rotor  will  then  begin  to  overshoot  the  vehicle  center  of 
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gravity.  It  has  been  taken  into  aeeount  that  this  is  a  traetor  type  tail  rotor  (starboard  side 
rotor).  The  traetor  type  rotor’s  thrust  is  in  the  opposite  direetion  of  the  vertieal  fin  and  the 
tail  rotor  inflow  stream  flows  towards  the  tail  fin. 

With  respeet  to  the  roll  rate  eontribution  to  the  hub  airspeed  in  the  y  direetion,  the 
absolute  value  of  the  roll  rate  is  neeessary.  As  the  vehiele  travels  in  the  positive  y 
direetion,  a  roll  rate  in  either  direetion  (port  or  starboard  roll)  will  eontribute  a  eomponent 
of  airspeed  in  the  negative  y  direetion  beeause  of  the  tail  rotor  position.  If  the  vehiele 
travels  in  the  negative  y  direetion,  the  roll  eomponent  will  have  the  effeet  of  inereasing 
the  tail  rotor  hub  airspeed  regardless  of  the  direetion  of  the  roll. 

The  normalized  rotor  inflow,  normal  to  the  tail  rotor  is  [1]; 


(3.36) 


The  tail  rotor  is  approximated  as  4.46  times  faster  than  the  main  rotor  speed 
(Q,^  =4.46Q^^)  and  is  ealeulated  as  =  778.2  rat/ /sec .  The  tail  rotor  hub  airspeed  in 
the  z  direetion  is  given  as; 


(3-37) 

where  is  the  vehiele  airspeed  in  the  z  direetion,  eompensated  for  wind  disturbanee, 
and  ( )  is  the  angular  airspeed  eomponent  due  to  a  pitehing  rate.  The  term  is  a 
wake  intensity  faetor  that  inereases  as  the  tail  rotor  beeomes  more  and  more  immersed  in 
the  main  rotor  wake.  The  rotor  wake  has  the  effeet  of  inereasing  the  airspeed  of  the  tail 
rotor  hub,  in  the  z  direetion.  The  term  is  the  main  rotor  indueed  veloeity.  In  low 
forward  airspeed  and  hovering  flight,  the  wake  intensity  faetor  and  grows  as  the 

tail  rotor  beeomes  more  and  more  immersed  in  the  main  rotor  wake  as  the  forward 
airspeed  inereases.  The  eondition  for  zero  wake  intensity  (near  hovering  flight)  is  given 
by[l]; 

(A,=0)  (3.38) 
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The  tail  rotor  is  fully  in  the  wake  when  [1]; 


V-  —w 

imr  a 


->gf  (^,=1.5) 


(3.39) 


in  whieh  ease,  K^-\.5  .  The  g,  and  g^  terms  are  determined  from  the  vehiele  geometry 
and  are  given  by  [1]; 


(3.40) 

(3.41) 

Equations  (3.38)  through  (3.41)  say  that  the  main  rotor  wake  does  not  impinge  on  the  tail 
rotor  airflow  until  the  vehiele  forward  airspeed  is  greater  than  6.25%  of  the  resultant 
main  rotor  inflow. 

Gavrilets  assumes  a  linear  growth  of  the  wake  intensity  with  inereasing  forward 
speed.  The  wake  intensity  is  given  by  [1]; 


Si  = 

Sf^ 


1  —  R  —  R 

Jr - mr - ^  =  .0625 


]  —  R  -\-  R 

jr - - ^  =  3.3125 


=l-5x 


Vjmr-W, 

Sf-g, 


■Si 


for 


V-  —w 

imr  a 


>gi 


and  Equation  (3.37)  ean  be  rewritten  as: 


V-  —w 

imr  a 


-.0625 


3.25 


The  tail  rotor  hub  airspeed  in  the  x  direetion  is  given  by: 

U,r  =  “a  +  14.7  +  hM 


(3.42) 


(3.43) 


(3.44) 
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Equation  (3.44)  states  that  the  tail  rotor  hub  airspeed  in  the  x  direetion  is  a  eombination 
of  the  vehiele  airspeed  u  and  the  piteh  and  yaw  rates.  Additionally,  it  states  that  the  piteh 
and  yaw  rates  eontribute  positive  airspeed  eomponents,  no  matter  the  sign  of  the  piteh 
and  yaw  rate.  For  example,  if  the  vehiele  is  traveling  in  the  positive  x  direetion,  a  piteh 
rate  in  either  direetion  (up  or  down)  will  eontribute  a  positive  airspeed  eomponent  to  the 
tail  rotor  hub,  in  the  x  direetion.  If  the  vehiele  travels  in  the  negative  x  direetion,  the 
positive  eontribution  of  the  piteh  or  yaw  rate  will  make  the  tail  rotor  airspeed  less 
negative  and  the  tail  rotor  begins  to  eateh  up  to  the  vehiele  eenter  of  gravity. 

With  all  of  the  tail  rotor  airspeed  eomponents  defined,  the  tail  rotor  hub  airspeed 
magnitude  is  then; 


The  tail  rotor  advanee  ratio  is  expressed  as  [1]; 


Mtr 


(3.45) 


(3.46) 


2,  Tail  Rotor  Inflow  and  Thrust 

To  determine  the  tail  rotor  inflow.  Equation  (3.11)  is  modified  to  suit  the  tail  rotor 
eonditions.  Equation  (3.11)  shows  that  the  main  rotor  inflow  is  the  sum  of  the  inflow  plus 
any  normal  airspeed  eomponent,  whieh  is  given  by  the  E^sin(a)  term.  In  the  tail  rotor 

ease,  the  tail  rotor  tip  path  plane  will  not  be  eanted  beeause  there  are  no  modeled  eyelies. 
Therefore,  the  airspeed  normal  to  the  tail  rotor  blades  will  be  given  simply  by  Equation 
(3.35).  Equation  (3.11)  now  beeomes: 

I  E 

Qi?  OR  (3.47) 
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Substituting  the  rotor  inflow  expression,  equation  (3.8),  into  (3.47)  gives  the  tail  rotor 
inflow  as: 


a 


The  Newton-Raphson  iteration  seheme  for  the  tail  rotor  is: 


fW 

f'W 


c. 


(3.48) 


(3.49) 


(3.50) 


(3.51) 


The  tail  rotor  inflow  is  determined  iteratively,  as  with  the  main  rotor,  using  Equations 
(3.49)  through  (3.51).  The  initial  estimate  (\)  is  the  normalized  tail  rotor  inflow  at 
hover.  The  tail  rotor  inflow  behaves  just  as  the  main  rotor  inflow  in  that  it  follows  the 
thrust  trend.  If  the  tail  rotor  thrust  inereases,  the  tail  rotor  inflow  inereases.  As  the  inflow 
inereases,  the  tail  rotor  thrust  will  settle  at  a  new  equilibrium  after  it  peaks.  This  is  more 
obvious  when  the  tail  rotor  thrust  is  expressed  as: 


T  - 


as 


^ped 


,3  2 


192.93/7 


^ped 


v3  2, 


(3.52) 


When  a  right  pedal  eommand  is  given  is  positive^,  using  the  previous  value 
of  the  inflow  (/I  at  hover)  will  eause  the  thrust  to  inerease.  The  new  value  of  the  tail 

rotor  is  fed  into  the  Netwon-Raphson  teehnique  and  gives  the  new  value  of  the  tail  rotor 
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inflow,  which  is  slightly  larger  than  the  previous  value.  The  tail  rotor  inflow  will  continue 
to  increase  until  the  rudder  input  has  settled.  Once  the  rudder  input  settles,  the  first  term 
in  (3.52)  settles  as  the  inflow  continues  to  increase.  This  causes  the  tail  rotor  thrust  to 
settle  after  it  attains  a  peak  value.  Figure  21  illustrates  the  tail  rotor  thrust  and  inflow 
behavior  to  a  step  in  tail  rotor  collective; 
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Figure  21 .  Tail  Rotor  Thrust  and  Inflow  Responses  to  Step  Collective 


The  tail  rotor  thrust  behaves  much  like  the  main  rotor  in  that  a  step  collective 
input  yields  the  eharacteristic  thrust  response  with  an  initial  rise,  followed  by  a  damping 
back  to  the  equilibrium  value. 

It  has  been  demonstrated  that  the  torque  produced  by  the  main  rotor  is 
- 6.297  Nm  .  Normalizing  this  by  the  tail  arm  ( 4,.  -  .9\m)  gives  the  force  that  the  tail 
rotor  needs  to  generate  to  counteract  the  main  rotor  torque.  This  gives: 
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‘-tr{hov) 


=  6.92  N 


Since  the  main  rotor  torque  produces  a  positive  yaw  (clockwise  rotation  of  the  vehicle  as 
viewed  from  above),  the  tail  rotor  force  at  hover  ( )  needs  to  point  in  the  positive  y 


direction  in  order  to  counteraet  the  main  rotor  torque. 

The  tail  rotor  eoefficient  of  thrust  can  be  determined  for  any  other  flight 
configuration  from  the  airspeed  and  angular  rate  conditions.  Once  the  tail  rotor  inflow  is 
determined,  the  tail  rotor  thrust  can  be  derived.  Padfield  gives  the  tail  rotor  thrust  as  [4,  p. 
141]; 


(3.53) 

where  the  coefficient  of  the  tail  rotor  thrust  is  [4,  p.  142]; 

- h - y  (3.54) 

Substituting  Equation  (3.54)  into  (3.53)  gives; 

(3.55) 

The  blockage  factor  accounts  for  thrust  losses  in  pusher  type  tail  rotors.  Since  this  is  a 
tractor  type  tail  rotor,  the  bloekage  factor  will  be  ignored  and  Equation  (3.55)  reduces  to; 

(3.56) 

The  net  tail  rotor  force  needed  to  counteract  the  main  rotor  torque  is  =  6.92  N  , 

which  is  the  same  as  the  thrust  required  and  gives  a  coefficient  of  thrust 
C*";)  =10.3954x10^1 
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At  hover,  the  normal  airflow  veloeity  is  equal  to  the  tail  rotor  inflow  veloeity  and 
Equation  (3.8)  beeomes  T  and  solving  for  the  indueed  veloeity  at  hover 

gives: 

^hcv  ^q29Amls 

i(tr) 

^hov 

pov  ^  ^  72.095  X 1 0^^ 


Solving  Equation  (3.52)  for  the  tail  rotor  oolleetive  required  at  hover  gives  Equation 
(3.57): 


^ped  - 


2C 

as 


T(tr)  2  _  //. 


ztr 

2 


3  2, 


(3.57) 


Solving  for  the  trim  tail  rotor  oolleetive  at  hover  gives 

as  2 
J"'”' =.19598  raJ  (11.23°) 

ped 

Tail  rotor  oolleetive  ranges  used  for  this  model  will  be  -15°  to  +25°  with  a  zero  oolleetive 
input  that  generates  the  trim  tail  rotor  oolleetive  of  1 1 .23°. 

The  tail  rotor  foroes  and  moments  are: 


Y,r=T,, 

Kr  =  Ur 


(3.58) 


3,  Yaw  Response  to  Tail  Rotor  Forces  and  Moments 

Eigure  22  illustrates  the  yaw  step  response  when  applying  the  tail  rotor  thrust  and 
moments  desoribed  in  Seotion  2. 
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Figure  22.  Tail  Rotor  Thrust  and  Inflow  Responses  to  Step  Collective 

The  yaw  angle  step  response  is  in  agreement  with  the  yaw  response  of  [1,  p.  1 12], 
as  expected.  The  yaw  angle  increases  in  a  counter-clockwise  direction  with  right  pedal 
input  and  damps  down  to  almost  zero  yaw  rate  after  the  pedal  input  is  taken  out.  The 
slight  drift  is  accounted  for  by  the  momentum  that  has  already  been  imparted  to  the  tail. 

D.  AERODYNAMIC  FORCES 

There  are  two  primary  aerodynamic  forces  acting  on  the  RW  UAV.  They  are  drag 
and  the  aerodynamic  forces  generated  by  the  main  rotor  airflow  around  the  UAV  body. 

1,  Determination  of  Drag  and  Main  Rotor  Tip-Path  Plane  Angle 

As  the  vehicle  approaches  forward  airspeeds  comparable  to  the  hover  induced 
velocity  (4  m/s),  drag  begins  to  increase  to  the  extent  that  it  has  to  be  compensated  for. 
Prouty  explains  that  the  tip  path  plane  angle  develops  as  a  direct  result  of  the  drag  forces 
that  act  on  the  vehicle.  For  a  vehicle  in  forward  flight,  the  drag  force  acts  in  the  negative 
X  direction.  In  order  to  balance  the  forces  in  the  x-direction,  a  longitudinal  cyclic 
command  must  be  put  in.  This  longitudinal  cyclic  input  causes  the  main  rotor  disk  to  tilt 
by  an  angle  (a)  with  respect  to  the  horizontal  plane.  As  the  disk  tilts,  the  thrust  vector 
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also  tilts  (thrust  vector  is  normal  to  the  disk)  and  creates  a  thrust  component  in  the 
positive  X  direction  that  will  counter  balance  the  drag  force. 

The  sum  total  of  all  drag  forces  on  the  vehicle  is  known  as  parasite  drag  and  is 
expressed  in  terms  of  an  equivalent  flat  plate  area,  /  .  According  to  Prouty  [3,  p.  132]; 

. .  .the  equivalent  flat  plate  area  is  the  frontal  area  of  a  flat  plate  with  a  drag 
coefficient  of  1,  which  has  the  same  drag  as  the  object  whose  drag  is  being 
estimated. 

Table  4  is  a  list  of  vehicle  components  for  Prouty’s  example  helicopter  along  with  initial 
estimates  for  the  T-Rex  Align.  The  baseline  for  estimates  of  the  T-Rex  Align  is  the 
fuselage  flat  plate  area.  Once  the  flat  plate  area  for  the  T-Rex  Align  is  estimated,  the 
conversion  factor  used  to  map  from  the  model  to  Prouty’s  example  is  determined  and 
used  for  all  the  other  components.  The  dimensions  used  for  estimating  the  vehicle 
fuselage  flat  plate  area  are  (.24x.24)  m,  this  gives  a  conversion  factor  of  9.93xl0^\ 
Components  that  contribute  very  little  parasite  drag  are  omitted 
[3,p.  132]: 


Component 

(Prouty’s  Example) 

Flat  Plate  Area  (ft^) 

X-Cell  60  Estimate  (m^) 

Fuselage 

5.8 

.05760 

Main  Rotor  Flub  and  Shaft 

7.0 

.06952 

Main  Landing  Gear 

1.2 

.01192 

Rotor-Fuselage  Interference 

1.3 

.01291 

Miscellaneous 

.5 

4.9655x10'^ 

Total 

15.8 

.1569 

Table  4.  Component-Wise  Flat  Plate  Area  Estimates 


The  total  parasite  drag  is  proportional  to  the  airspeed  and  is  expressed  as  [3,  p.  132]; 

(3.59) 

The  main  rotor  tip-path  plane  is  approximated  as  [3,  p.  133]; 

a  =  -573  ^  (deg)  (3.60) 

GW. 
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This  angle  is  the  main  rotor  tip  path  plane  angle,  with  respeet  to  the  horizontal 
plane,  that  is  required  to  keep  the  vehiele  in  trim  flight,  for  a  given  forward  veloeity.  It 
will  be  shown  later,  how  the  fuselage  dynamies  affeet  this  angle  as  well  as  how  it  relates 
to  the  pilot  eyelie  input.  Sinee  the  thrust  veetor  is  normal  to  the  tip-path  plane,  it  follows 
that  the  thrust  in  the  z-direetion,  ,  is: 


^MR^Tcos[a)  {N)  (3.61) 

The  foree  in  the  x-direetion  due  to  the  main  rotor  is: 

XMR^Tsm{a)  {N)  (3.62) 

Balaneing  forees  in  the  x  direetion  gives: 

Tsm{a)-f^V^  (3.63) 

Solving  for  the  trim  thrust  required  in  this  flight  eonfiguration  gives: 

{Forward  Flight :  k  >  v,) 

fEyl 

r  =  (iV)  (3.64) 

sin(^aj 

2 

^  _ ^ihover _ 

+  2Tv,.  sin  a  +  vf 

Again,  this  is  the  trim  eondition  thrust  required  to  balanee  the  main  rotor  forees,  the 
vehiele  weight,  and  the  parasite  drag. 

2.  Fuselage  Forces 

According  to  Gavrilets,  the  rotor  downwash  is  deflected  by  the  forward  and  side 
velocity,  when  near  the  hover  flight  regime.  The  deflection  of  the  rotor  downwash  creates 
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a  force  that  acts  opposite  the  direction  of  movement.  The  fuselage  drag  forces  in  the  X 
and  Y  directions  are  given  by  [1,  p.  47] 

(3.65) 

>7,.=5r|p'’,...v  (3.66) 

where  5'/“  and  5'f'^  are  projected  cross  sectional  areas  of  the  fuselage  on  the  YZ  and  XZ 
planes. 

In  forward  flight  (when  the  translational  velocity  is  greater  than  the  main  rotor 
inflow  velocity),  the  fuselage  drag  forces  are  expressed  as  [3] 


^fus-S^^pU^U 

(3.67) 

(3.68) 

Where  Ug  is  the  trim  airspeed.  Gavrilets  [1]  gives  the  tail  rotor  hub  velocity  magnitude 

as: 

+  '^ihovf 

(3.69) 

and  Equations  (3.67)  and  (3.68)  can  be  expressed  as  [1]: 

^fus  —\pS{"\V^ 

(3.70) 

Yfus--\pS^‘‘\V^ 

(3.71) 

The  fuselage  drag  force  in  the  Z  direction  is  expressed  as  [1]; 

(3.72) 


Mettler  assumes  that  «  2.2S{'‘\  « 1.5^/“ ,  where  is  determined  by  the 

vehicle  cross  sectional  area  with  respect  to  the  YZ  axis,  when  the  vehicle  was  operating 


48 


at  15.4  m/s  forward  air  speed.  This  resulted  in  5'/“  =  .l  m^.  This  initial  estimate  will  also 
be  used  for  modeling  T-Rex  Align  UAV.  The  remaining  cross  sectional  areas  are; 

5'f  *  =  .22  rn 
=  .15 

In  summary,  the  main  rotor  thrust  is  calculated  by  using  Equation  (3.1),  which 
relates  the  operator  commanded  collective  (6*^),  the  vehicle  velocity  states,  the  main 

rotor  aerodynamic  and  dimensional  parameters,  and  the  main  rotor  inflow  (T.).  The 

Newton-Raphson  iteration  technique  is  used  to  compute  the  convergent  inflow  at  a  rate  of 
100  Hz.  The  main  rotor  flapping  dynamics  is  modeled  as  a  second  order,  linear  system 
whose  inputs  are  commanded  flap  angles,  the  pitch  rate,  and  roll  rate.  The  tail  rotor  was 
modeled  like  the  main  rotor,  with  the  exception  that  there  are  no  tail  rotor  cyclics.  This 
ensures  that  the  tail  rotor  thrust  always  acts  normal  to  the  helicopter  fin  plane.  Also,  the 
tail  rotor  has  its  own  tail  rotor  hub  velocities,  which  are  not  necessarily  the  same  as  the 
center  of  gravity  velocity  states.  The  forces  and  moments  on  the  vehicle  center  of  gravity 
are  computed  using  Equations  (3.33)  and  (3.34).  These  forces  and  moments  are  summed 
to  provide  force  and  moment  components  to  the  6-DOF  block  set. 

The  following  chapter  will  go  through  the  linearization  of  the  nonlinear  model  in 
order  to  uncover  the  state  boundaries  for  which  the  controller  design  will  be  useful. 
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IV.  LINEARIZATION  OF  THE  NONLINEAR  ROTARY-WING 

MODEL 


The  6-DOF  RW  UAV  model  is  highly  nonlinear.  Linearizing  the  model  about  a 
eertain  operating  point  makes  the  eontrol  design  process  simpler.  The  goal  is  to  obtain  the 
necessary  linear  models  that  describe  the  helicopter  dynamics  and  to  design  the  PID 
controllers  that  will  minimize  the  flight  trajectory  errors  during  SIL  testing. 

Matlab  has  several  linearization  tools  that  can  be  applied  to  most  nonlinear 
models.  The  tools  that  were  applied  for  this  linearization  task  are  the  functions  “linmod” 
and  “n4sid.”  The  linmod  function  generates  a  state  space  model,  given  the  model  name, 
initial  state  vector,  and  inputs.  The  “n4sid”  function  is  used  to  obtain  linear  models  for 
dynamics  that  need  to  be  identified  through  system  identification.  Please  see  Appendices 
C  and  D  for  Matlab  script  implementation  of  these  functions. 

A.  TRIMMING  OF  THE  NONLINEAR  MODEL 

The  trim  solution  is  given  by  the  control  inputs 

W.  =K„  <  Kol  KeA 

required  to  maintain  the  helicopter  at  specified  states 

x^=\u  V  w  (/)  6  y/  p  q 

Trimmed  flight  requires  that  the  rate  of  change  (of  magnitude)  of  the  aircraft’s 
state  vector, 

Xe  =0 

and  the  sum  of  the  forces  and  moments  on  the  aircraft  are  zero.  Once  an  operation  point 
or  steady  states  are  chosen,  the  aircraft  will  be  trimmed.  That  is,  the  trim  inputs  required 
to  maintain  the  aircraft  in  the  specified  operating  condition  will  be  determined.  The  trim 
inputs  and  states  are  determined  because  this  allows  the  use  of  small  perturbation  theory 
to  make  estimates  of  the  helicopter  states  about  a  trim  operating  point.  Small  perturbation 
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theory  says  that,  during  perturbed  motion,  the  helieopter  behavior  ean  be  deseribed  as  a 
perturbation  from  the  trim,  written  as  [4,  Pg.  208]; 

x  =  x^+5x 

The  operation  point  that  was  ehosen  for  trimming  is  the  hovering  flight  eondition 
where  all  states 

x^=\u  V  w  (/)  6  y/  p  q 

are  zero.  Appendix  D  outlines  the  eommands  given  to  find  the  trim  settings  for  any 
speeified  set  of  state  eonditions.  The  desired  states  are  speeified  and  the  error  toleranees 
are  defined.  At  first,  the  nonlinear  model  response  is  simulated  with  an  initial  estimate  of 
the  trim  inputs  required.  The  state  values  at  the  end  of  the  simulation  are  then  taken  and 
eompared  to  the  desired  values  and  the  errors  are  eomputed.  The  eomputed  errors  are 
then  used  to  adjust  the  appropriate  eontrol  input  and  the  simulation  is  run  again.  This 
proeess  is  eondueted  iteratively  until  the  adjusted  trim  eonditions  ean  maintain  the 
desired  helieopter  states  within  the  defined  toleranees. 

B,  LINEARIZATION  OF  THE  LATERAL  AND  LONGITUDINAL 
DYNAMICS 

The  lateral  and  longitudinal  dynamies  are  generated  by  the  main  rotor  flapping 
response  along  the  y  and  x  axes,  respeetively.  A  stiek  eommand  (left-right)  for  a  lateral 
flap  angle  ,  for  example,  will  eause  the  lateral  eyelie  servo  to  rotate  and  generate  a 

main  rotor  flap  angle  eommand.  This  tilting  of  the  rotor  in  the  y  direetion  generates  a 
foree  in  the  y  direetion  as  well  as  a  rolling  moment,  whieh  eauses  the  UAV  to  roll.  A 
stiek  eommand  (forward-aft)  for  a  longitudinal  flap  angle  ,  on  the  other  hand,  will 

eause  the  longitudinal  eyelie  servo  to  rotate  and  generate  a  main  rotor  flap  angle 
eommand  in  the  longitudinal  direetion  (+/-x  direetion).  The  tilting  of  the  rotor  in  the  x 
direetion  generates  a  foree  in  the  x  direetion  as  well  as  a  pitehing  moment. 

Onee  the  hover  trim  eonditions  have  been  determined,  as  deseribed  in  seetion  A, 
the  lateral  and  longitudinal  dynamies  linear  models  will  be  obtained  by  using  the 
Simulink  funetion  “linmod.”  The  linmod  funetion  extraets  the  state  spaee  matriees  A,  B, 

C,  and  D  given  the  equilibrium  state  and  input  veetors,  and  ,  respeetively.  Onee  the 
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state  space  representation  is  derived,  the  minimum  realization,  single  input  single  output 
(SISO)  model  is  extracted.  The  minimum  realization  reduces  the  model  order  and 
complexity  for  states  that  contribute  little  to  no  dynamic  behavior  for  the  degree  of 
freedom  in  question.  The  details  on  the  commands  used  to  generate  the  lateral  and 
longitudinal  models  are  in  Appendices  E  and  F. 

Figure  23  illustrates  the  6-DOF  nonlinear  model  with  four  command  inputs  and 
nine  state  outputs.  Figure  24  illustrates  more  of  the  model’s  internal  structure. 


T-Rex  Align  600  Model 
Ver.  1.1 

Based  on  original  work  created  by 
Capt  Ricardo  E  .  Miranda  for  the  Naval  Postgraduate  School 

Monterey  ,  CA 
2009 
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Figure  23.  Nonlinear  6-DOF  T-Rex  Align  600  Model 
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Figure  24.  Nonlinear  6-DOF  T-Rex  Align  Model  (Internal  Strueture) 

The  file  “RunModel.m”  (Appendix  D)  is  run  in  order  to  generate  step  input 
responses  and  to  generate  linear  models  for  all  six  degrees  of  freedom.  Onee  all  of  the 
linear  models  are  obtained,  the  nonlinear  step  responses  are  eompared  to  the  linear  model 
step  responses.  A  good  linear  model  provides  similar  responses  to  the  nonlinear  model 
for  a  short  period  of  time,  before  the  linear  response  begins  to  diverge  from  the  nonlinear 
response.  Appendix  G  illustrates  the  internal  structure  of  the  T-Rex  Align  600  linear 
model. 

Figure  25  illustrates  the  open-loop  step  response  of  the  nonlinear  model  to  a  2° 
longitudinal  step  command  in  flap  angle  as  well  as  the  linear  model’s  response. 
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Figure  25.  Longitudinal  Dynamics  /  Linear  Model  Validation 

Figure  25  illustrates  that  the  linear  model  represents  the  nonlinear  longitudinal 
model  quite  accurately  in  the  velocity  u.  The  pitch  response  of  the  helicopter  is  not 
directly  controllable  and  is  designed  with  a  high  degree  of  dynamic  stability.  Controller 
design  is  designed  around  the  longitudinal  velocity  model  obtained.  Figure  25  also 
illustrates  that  the  linear  model  is  accurate  in  representing  the  surge  velocity  of  the 
vehicle  up  through  1.5  m/s.  This  indicates  that  the  helicopter  model  can  be  disturbed  from 
hover,  out  to  1.5  m/s,  without  departing  from  its  linear  region.  This  observation  makes  it 
possible  to  design  a  flight  path  controller  near  the  hover  flight  regime.  The  input  to  output 
transfer  functions  and  eigenvalues  for  the  principal  dynamics  are  listed  in  Appendix  FI. 

Figure  26  illustrates  the  step  response  of  the  nonlinear  model  to  a  lateral  step 
command  of  2°  in  flap  angle.  Also  illustrated  are  the  linear  model  responses. 
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Lateral  Step  Input  (deg) 
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Figure  26.  Lateral  Dynamies/L inear  Model  Validation 


The  vehicle  behaves  differently  when  disturbed  in  the  lateral  direction,  as 
compared  to  the  longitudinal  responses.  This  is  because  the  helicopter  has  a  tendency  to 
rotate  in  the  direction  of  lateral  disturbance,  in  the  absence  of  tail  rotor  thrust 
compensation.  As  the  vehicle’s  yaw  angle  increases  beyond  90°,  the  velocity  in  the  y 
direction  begins  to  decrease,  as  illustrated  in  the  second  figure  of  Figure  26.  This  occurs 
around  4.5  seconds  simulation  time.  This  explains  the  divergence  of  the  nonlinear 
response  in  v  to  the  linear  model  response.  This  indicates  that  the  linear  model  may  be 
valid  well  through  a  disturbance  of  1  m/s  from  the  hover  state.  The  input  to  output 
transfer  functions  for  the  principal  dynamics  are  listed  in  Appendix  FI. 
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c. 


LINEARIZATION  OF  THE  HEAVE  DYNAMICS 


Attempts  to  linearize  the  heave  dynamies  using  the  same  teehniques  outlined  in 
seetion  B  did  not  provide  aeeurate  linear  models.  While  the  lateral  and  longitudinal 
dynamies  are  governed  by  the  main  rotor  flapping,  the  heave  dynamies  are  governed  by 
the  nonlinear  thrust  response  to  collective  commands.  The  linearization  approach  that 
was  used  to  model  the  heave  response  is  system  identification.  The  goal  is  to  obtain  a 
system  identification  model  by  providing  the  “n4sid”  function  input  and  output  data.  It 
should  also  be  mentioned  here  that  the  heave  dynamics  in  the  climb  and  descent  are  quite 
different,  as  explained  in  the  modeling  sections  earlier.  Therefore,  two  heave  models 
were  obtained,  one  for  climb  and  one  for  descent. 

Figure  27  illustrates  the  climbing  heave  step  response  for  the  linear  and  nonlinear 
model  responses. 


MR  Step  Input  (deg) 


Velocity  (m/s) 


Figure  27.  Fieave  Dynamics/Linear  Climb  Model  Validation 
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Figure  27  illustrates  the  linear  elimbing  heave  model  is  accurate  in  the  principal 
heave  dynamics  and  the  that  the  linear  model  is  accurate  in  representing  the  nonlinear 
model  up  to  3  m/s  climbing  speed,  given  the  2°  step  collective  input.  In  the  flight  path 
controller  design,  collective  commands  will  only  be  sufficient  to  maintain  the  vehicle  at  a 
constant  altitude,  and  therefore  heave  velocities  will  be  well  within  the  3  m/s. 

Figure  28  illustrates  the  linear  descent  heave  model  response  compared  to  the 
nonlinear  heave  descent  response. 


MR  Step  Input  (deg) 
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Figure  28.  Fleave  Dynamics/Linear  Descent  Model  Validation 

Figure  28  illustrates  that  the  linear  heave  descent  model  obtained  diverges  slightly 
from  the  nonlinear  model  at  around  3  m/s  descent.  Inside  of  this  boundary,  the  linear 
model  represents  the  nonlinear  dynamics  very  closely.  The  heave  transfer  functions  for 
climb  and  descent  are  listed  in  Appendix  H. 
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D.  LINEARIZATION  OF  THE  YAW  DYNAMICS 


The  yaw  dynamics  are  generated  by  the  combination  of  the  main  rotor  torque  on 
the  UAV  body  and  the  tail  rotor  thrust  compensation.  A  stick  command  (left-right)  for  a 
tail  rotor  cyclic  command  will  cause  the  tail  rotor  thrust  to  increase  or  decrease, 
depending  on  the  type  of  command.  The  changing  of  the  tail  rotor  thrust  generates  a  force 
on  the  UAV  tail,  which  causes  the  UAV  to  yaw.  The  linearization  process  that  was 
applied  for  the  lateral  and  longitudinal  dynamics  cannot  be  applied  for  the  yaw  degree  of 
freedom.  The  approach  that  was  taken  here  is  system  identification  using  a  data  driven 
modeling  approach.  The  goal  was  to  use  simulation  data  from  the  nonlinear  model  to 
obtain  a  linear  model  of  the  yaw  dynamics  from  the  commanded  tail  rotor  collective  to 
the  yaw  angle  y/ .  Figure  29  illustrates  the  nonlinear  yaw  response  and  the  linear  model 
yaw  response. 
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Figure  29.  Yaw  Dynamics  Model  Validation 
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Figure  29  illustrates  that  the  linear  yaw  dynamics  model  obtained  is  accurate  in 
the  principal  yaw  dynamics,  y/ .  The  yaw  dynamics  transfer  function  is  listed  in 
Appendix  H. 

In  summary,  the  nonlinear  model  was  linearized  about  the  longitudinal  and  lateral 
modes  using  the  Matlab  function  “linmod,”  while  the  heave  and  yaw  modes  were 
linearized  using  system  identification;  the  system  identification  function  “n4sid” 
facilitated  this  process.  For  a  2°  flap  angle  step  command,  the  linear  model  is  accurate 
when  compared  to  the  nonlinear  responses  in  the  longitudinal  and  lateral  modes.  The  step 
responses  were  accurate  up  to  1  m/s  and  they  also  indicate  that  linear  models  will  provide 
accuracy  well  through  1  m/s.  This  was  not  verified  and  for  this  reason,  flight  path 
tracking  performance  of  the  PD  gains  to  be  derived  will  remain  within  this  flight  regime. 
The  linear  model  for  heave  is  also  accurate  up  to  3  m/s,  in  both  climb  and  descent.  The 
linear  model  for  the  yaw  mode  is  also  very  accurate  for  a  1°  step  command  in  tail  rotor 
collective. 

The  next  chapter  will  focus  on  the  development  of  the  closed  loop  system  and 
tuning  of  the  PD  controllers  for  the  longitudinal,  lateral,  heave,  and  yaw  modes. 
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V.  CONTROL  DESIGN 


A,  CLOSED  LOOP  SYSTEM  STRUCTURE 

Figure  30  illustrates  the  nonlinear  model  with  a  lead  Proportional  Derivative  (PD) 
eompensator.  The  PD  eompensator  reeeives  the  error  signal,  whieh  is  the  differenee 
between  the  referenee  signal  (desired  position)  and  unit  feedbaek  output  (aetual  position). 


Figure  30.  Closed-Loop  System  with  Lead  PD  Compensation 

The  primary  degrees  of  freedom  to  be  eontrolled  are  longitudinal,  lateral,  heave, 
and  yaw.  This  is  aehieved  by  implementing  four  Proportional  Derivative  (PD)  eontrollers 
for  the  aforementioned  modes. 

B,  DETERMINATION  OF  PD  CONTROLLER  GAINS 

The  Ziegler-Niehols  tuning  method  outlined  in  [8,  p.  673]  was  applied  in  order  to 
obtain  appropriate  PD  gains.  The  Ziegler-Niehols  tuning  method  starts  with  all  gains  at 
zero  and  the  test  gain  is  inereased  until  sustained  oseillations  are  obtained.  A 
sustained  oseillation  is  eharaeterized  by  an  oseillatory  output  that  is  eonstant  in  amplitude 
and  period  [8,  p.  672].  Onee  the  oseillation  period  and  test  gain  this  point  are 
known.  Table  5  is  applied  in  determining  all  of  the  initial  estimate  gains. 
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Control  Type 

Kp 

K, 

Ko 

P 

.5Kp 

PI 

A5Kp 

llKpIP^ 

PID 

.6Kp 

2Kp/P^ 

KpPJS 

Table  5.  Ziegler-Nichols  PID  Gain  Estimates  (After  [8]) 

To  determine  the  initial  longitudinal  PD  gains,  a  step  input  (input  port  1)  is 

provided  to  the  closed  loop  system  illustrated  in  Figure  30.  With  no  longitudinal  PD 
compensation  {Kp  =  =  0)  ,  a  reference  step  input  in  the  longitudinal  gives  no  output 

position  signal  x  (output  port  1),  since  the  error  signal  into  the  channel  is  zero.  The  test 
gain  Kp  is  then  increased  to  Kp=  A  and  the  output  position  analyzed. 

Figure  3 1  illustrates  the  position  response  to  a  reference  position  step  input  with 
Kp=A  and  oscillation  period  =11  sec.  Figure  31  illustrates  that  a  sustained 

oscillation  has  been  obtained  with  Kp  = .  1 . 


Desired  X  Position  (m) 


Figure  3 1 .  Position  Response  {Kp=A)l  PD  Tuning 
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Figure  3 1  illustrates  the  error  signal  along  with  the  control  input  u  that  is 
generated  from  the  proportional  gain.  This  shows  that  the  position  response  is  oscillating 
about  the  desired  response  of  1  meter  at  the  same  frequency  as  the  control  input  u.  It  also 
shows  that  the  error  signal,  which  drives  the  control  input,  is  too  large  and  is  why 
insufficient  control  effort  is  being  applied  in  this  case. 

Using  the  PID  row  of  Table  5,  without  integral  control,  gives  the  initial 
longitudinal  PD  gain  estimates  of  Kp  =  .06  and  =  .825  .  An  integral  gain  is  not  used 

since  doing  so  made  the  closed  loop  response  go  unstable.  Figure  32  illustrates  the 
position  response  with  the  initial  PD  gain  estimates. 
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Figure  32.  Position  Response  {Kp=  .06  =  .825 )/  Longitudinal  Mode 

Figure  32  shows  that  the  error  starts  off  as  1  meter  and  damps  down  exponentially 
towards  zero.  This  is  generated  by  an  initial  control  impulse  at  1  second  with  a  peak  of  .8. 
The  total  time  over  which  the  control  input  is  applied  is  roughly  .35  seconds. 

Figure  32  also  shows  that  the  response  to  a  step  input  of  1  meter  is  too  slow.  The 
desired  response  should  yield  a  settling  time  inside  of  3  seconds  with  an  overshoot  of  less 
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than  5%.  Tuning  the  gains  further  gives  Kp  =3.0 and  =1.6  with  the  step  responses 

illustrated  in  Figure  33.  The  lateral  dynamies  are  very  similar  to  the  longitudinal  and 
using  the  same  gains  provided  similar  results. 


Desired  X  Position  (m) 


Figure  33.  Position  Response  ( =  3.0  =  1 .6 )/  Longitudinal  Mode 


Figure  33  illustrates  that  the  error  between  the  desired  and  aetual  position  damps 
down  to  within  2%  of  zero  error  1.11  seeonds  after  the  referenee  signal  is  eommanded. 
The  determined  PD  gains  yield  a  settling  time  less  than  3  seeonds  and  an  overshoot  less 
than  3%.  Figure  34  shows  the  veloeity  response  along  with  the  position. 
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X  Position  Response  (m) 


time  (sec) 

Figure  34.  Velocity  Response  ( =  3.0  =  1 .6 )/  Longitudinal  Mode 


The  same  tuning  procedure  is  performed  for  the  heave  PD  controller  gains.  Figure 
35  illustrates  the  Z  velocity  response  along  with  the  PD  compensator  error  signal  and 
resultant  control  signal  u  . 
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Figure  35.  Position  Response  {Kp=  .45  =  .2 )/  Heave  Mode 
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The  heave  response  gives  a  settling  time  less  than  2  seeonds  with  no  overshoot. 

The  yaw  mode  controller  gains  yield  the  responses  illustrated  in  Figure  36. 


Tail  Rotor  Step  Collective  (deg) 


Figure  36.  Yaw  Angle  Response  {Kp=.\  =  .05 )/  Yaw  Mode 


The  yaw  angle  response  gives  excessive  overshoot  but  this  can  be  compensated 
for  in  the  future  by  including  a  yaw  rate  feedback.  For  now,  this  will  be  sufficient  to 
provide  a  steady  orientation.  Table  6  outlines  the  PD  gains  for  all  four  modes. 


MODE 

Kp 

Ko 

Longitudinal 

3.0 

1.6 

Lateral 

3.0 

1.6 

Fleave 

.45 

.2 

Yaw 

1.1 

.25 

Table  6.  PD  Gains  at  Flover 
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C.  FLIGHT  PATH  TRACKING  PERFORMANCE  /  SQUARE  PATH 

The  PD  controller  gains  listed  in  Table  6  are  tested  on  the  nonlinear  model  by 
feeding  the  closed  loop  system  a  table  of  reference  inputs  that  represent  position 
commands  in  the  three  dimensions.  The  following  will  analyze  the  controller  tracking 
performance. 

1,  Position  Tracking  Performance 

The  flight  path  trajectory  that  was  set  up  for  this  test  is  a  4mx4m  square  at  a  2- 
meter  altitude.  The  initial  position  of  the  model  is  set  at  the  origin  and  the  first  leg  of  the 
square  path  is  along  the  y  -axis.  Figure  37  illustrates  the  x  and  y  position  response  of  the 
model  versus  the  reference  path. 
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Figure  37.  Flight  Path  Tracking  Performance/Square  Path 
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Figure  37  illustrates  that  the  longitudinal  and  lateral  mode  PD  controller  performs 
well.  The  maximum  position  error  observed  is  3.1  cm  and  is  more  easily  identified  in 
Figure  38,  which  is  a  zoomed  in  version  of  Figure  37,  for  the  first  leg  of  the  track. 


Trajectory  Tracking  Performance 


Figure  38.  Flight  Path  Tracking  Performance/Error 


Figure  39  illustrates  the  nonlinear  model’s  flight  path  in  three  dimensions.  Figure 
40  illustrates  the  flight  path  trajectory  with  increased  vertical  position  resolution,  which 
shows  that  the  maximum  vertical  position  error  is  about  1 .4  cm. 
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Figure  39.  3D  Flight  Path  Tracking  Performance/Square  Path 
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Square  Trajectory  Tracking  Performance 


Figure  40.  3D  Flight  Path  Tracking  PerformanceA^ertical  Error 


2,  Yaw  Control  Performance 

Figure  41  illustrates  the  yaw  response  (in  degrees)  for  a  commanded  yaw  angle  of 
zero  degrees.  Initially,  as  the  helicopter  is  climbing  to  2  meters,  the  tail  boom  is  rotating 
clockwise  because  of  the  increasing  main  rotor  moment  on  the  body  and  explains  the 
increasing  yaw  angle  for  the  first  three  to  four  seconds.  The  remaining  25  seconds  of  the 
response  illustrated  in  Figure  41  represents  the  first  leg  where  the  helicopter  model  is 
traveling  from  left  to  right,  along  the  inertial  frame’s  positive  y  -axis.  The  reference 
positions  are  commanded  in  a  way  that  causes  the  helicopter  to  surge  in  1  meter 
increments.  Then,  the  helicopter  slows  down  as  it  approaches  its  final  position,  and 
hovers  until  the  next  1  meter  command  is  given.  This  explains  the  4  yaw  angle  spikes 
illustrated  in  Figure  41.  As  the  model  travels  in  its  sideways  trajectory  (1  meter  at  a  time), 
the  tail  boom  has  the  tendency  to  rotate  clockwise  and  explains  the  increasing  yaw  angle. 
As  the  yaw  angle  increases,  the  PD  controller  applies  control  effort  and  brings  the  yaw 
angle  back  to  zero,  before  the  model  begins  another  1  meter  surge  in  the  positive  y 
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direction.  The  yaw  angle  spikes  can  be  reduced,  by  applying  a  yaw  rate  feedback  to  the 
nonlinear  model.  Figure  42  illustrates  the  yaw  angle  responses  for  the  entire  flight  path 
tracking  simulation  and  shows  that  the  yaw  angle  errors  are  smaller,  for  the  second  and 
fourth  legs  of  the  flight  path,  because  of  the  flight  path  orientation  with  respect  to  the 
helicopter  tail  boom. 


Yaw  Angle  Response 


Figure  41.  Yaw  Angle  Response/First  Leg 


Yaw  Angle  Response 


Figure  42.  Yaw  Angle  Response/Entire  Path 
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3,  Euler  Angle  Responses 

This  model  does  not  implement  a  direct  control  scheme  for  the  UAV  orientation; 
rather,  it  relies  on  the  inherent  stability  that  has  been  modeled.  The  UAV  has  been 
modeled  so  that  an  angle  disturbance  will  be  compensated  for  by  the  hub  spring  force  as 
much  as  possible  without  having  undesired  effects.  Figure  43  illustrates  the  model 
position  response  in  the  y  direction  along  with  the  roll  angle  response,  for  the  first  leg  of 
the  simulation. 
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Figure  43.  Roll  Angle  Response/First  Leg 


Figure  43  illustrates  that  with  each  successive  surge  in  the  y  direction,  the  roll 
angle  increases  in  steps.  With  each  surge,  the  roll  angle  increases  then  damps  back  to  a 
smaller  roll  angle  that  is  larger  than  the  initial  roll  angle.  Also  shown  is  the  pitch  angle 
which  is  only  slightly  increasing  (positive  pitch  means  nose  up)  due  to  rotor  flap  angle 
coupling.  Figure  44  illustrates  the  y  position  response  and  the  roll  angle  responses  for  the 
entire  simulation.  It  is  evident  that  the  roll  angle  begins  to  decrease  back  towards  zero  as 
the  UAV  travels  in  the  -y  direction. 
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Y  Position  Response  (m) 
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Figure  44.  Roll  Angle  Response/Entire  Simulation 

Figure  45  shows  similar  results  for  the  piteh  angle  responses.  As  the  x  position 
increases,  the  pitch  angle  becomes  more  negative  (nose  down).  The  pitch  angle  begins  to 
recover  as  the  UAV  flies  backward  along  the  x  -axis. 


X  Position  Response  (m) 


Pitch  Angle  Response  (deg) 


Figure  45.  Pitch  Angle  Response/Entire  Simulation 
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D.  FLIGHT  PATH  TRACKING  PERFORMANCE  /  FIGURE  9 

A  figure-8  reference  flight  path  is  generated  with  user  specified  radius,  number  of 
rotations  for  each  loop,  and  desired  flight  time  for  each  complete  circle.  Initially, 
the  model  begins  at  ground  level  and  climbs  to  5  meters  and  once  positioned  at 
[x  y  z]=[2  0  5],  the  figure-8  path  tracking  routine  begins.  The  radius  chosen  for  this 
demonstration  is  1  meter  with  the  initial  position  of  [2  0  5].  The  desired  time  to  traverse 
the  360°  for  each  loop  is  defined  as  10  seconds. 

1,  Position  Tracking  Performance 

Figure  46  illustrates  the  top  view  of  the  helicopter  figure-8  trajectory  versus  the 
commanded  reference  trajectory. 


Trajectory  Tracking  Performance 


Figure  46.  Flight  Path  Tracking  Performance/Figure  9  (Top  View) 

Using  the  same  PD  gains  derived  earlier  gives  acceptable  tracking  performance 
with  a  maximum  position  error  of  3.6  cm.  Figure  47  illustrates  the  three  dimensional 
flight  trajectory  with  the  initial  climb  command. 
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Trajectory  Tracking  Performance 


Figure  47.  Flight  Path  Tracking  Performance/Figure  9  (3D) 


Figure  48  is  a  zoomed-in  version  to  show  more  resolution  in  the  vertical  axis.  The 
maximum  error  in  the  vertical  position  is  about  1 .7  mm  after  an  initial  maximum  error  of 
8.6  mm. 


Vertical  Error 


Figure  48.  Flight  Path  Tracking  Performance/Figure  9  (Vertical  Position  Error) 
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Figure  49  illustrates  the  velocity  responses  u  and  v. 


Velocity  Responses  (m/s) 


Figure  49.  Flight  Path  Tracking  Performance/Figure  9  (Velocity  Responses) 

2,  Yaw  Control  Performance 

Figure  50  illustrates  the  yaw  angle  response  to  a  commanded  angle  of  zero 
degrees.  The  swaying  of  the  tail  boom  is  natural  in  that  the  tail  has  the  tendency  to  rotate 
anytime  the  vehicle  turns  in  one  direction  or  another.  If  the  vehicle  turns  right,  the  tail 
will  rotate  clockwise  and  if  the  vehicle  turns  left,  the  tail  will  rotate  counter  clockwise.  In 
this  ligure-8  trajectory,  the  vehicle  is  constantly  turning  and  is  why  the  tail  is  constantly 
oscillating.  In  this  simulation,  the  tail  does  not  sway  more  than  4.5°  once  established  on 
the  Figure-8  trajectory.  As  mentioned  earlier,  however,  the  tail  boom  yaw  angle  error  can 
be  reduced  significantly  by  implementing  a  yaw  rate  feedback  in  future  control  projects. 


75 


Yaw  Angle  Response  (deg) 


F igure  50.  Yaw  Angle  Response/F igure  9 


3,  Euler  Angle  Responses 

Figure  5 1  illustrates  the  roll  angle  response  with  respeet  to  the  y  position. 


Y  Position  Response  (m) 
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Figure  5 1 .  Roll  Angle  Response/Figure  9 
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From  about  12  seconds  to  15  seconds,  Figure  51  illustrates  the  vehicle  y  position. 
In  this  time  frame,  the  vehicle  starts  at  y  =  lm  and  flies  towards  y  =  0.  Analyzing  the 
bottom  figure  of  Figure  51,  along  the  same  time  frame,  shows  that  the  roll  angle  is 
decreasing  (rolling  towards  the  left  side)  and  is  consistent  with  the  desired  roll  angle 
response.  In  traveling  from  y  =  -1  mback  to  y  =  0,  the  roll  angle  is  increasing  (rolling  to 
the  right).  These  observations  are  consistent  with  the  expected  roll  angle  responses. 

Figure  52  illustrates  the  pitch  angle  response  along  with  the  corresponding  x 
position. 


X  Position  Response  (m) 
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Figure  52.  Pitch  Angle  Response/Figure  9 

Figure  52  illustrates  that  as  the  vehicle  position  along  the  x  axis  decreases 
(backwards  flight),  the  pitch  angle  increases  (increasingly  nose  up  attitude).  Conversely, 
as  the  vehicle  travels  along  the  +x  -  axis  (forward  flight),  the  nose  pitches  down  for  a 
decreasing  pitch  angle. 
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E, 


CONCLUSIONS  ON  TRACKING  PERFORMANCE 


The  position  tracking  performance  has  given  nice  results  with  maximum  errors  of 
3.6  cm  for  translational  motion  and  8  mm  for  vertical  position.  The  yaw  angle  errors  are 
expected  given  that  no  yaw  rate  feedback  is  modeled  at  this  point.  The  actual  plant  does 
have  a  yaw  gyro  that  provides  rate  feedback  to  the  sensors.  Mettler’s  work  [2]  is  a  good 
source  for  obtaining  yaw  rate  feedback  models.  The  angle  responses  are  generally  correct 
in  that  main  rotor  flap  commands  generate  the  roll  and  pitch  responses  that  are  expected. 
Flight  in  the  ±x  direction,  for  example,  generates  the  anticipated  pitch  angle  responses. 
Conversely,  flight  in  the  ±y  direction  generates  the  anticipated  roll  angle  responses. 
Getting  the  specific  quantitative  RW  UAV  roll  and  pitch  angle  responses  will  require 
system  identification  of  the  MAE  Department  miniature  helicopter  so  that  the  flapping 
spring  derivatives  and  can  be  determined.  As  explained  earlier,  these  derivatives 

determine  the  way  that  aerodynamic  moments  from  the  main  rotor  are  cancelled  by 
restoring  moments  generated  by  the  blade  centrifugal  forces  and  by  damping  moments 
provided  by  the  main  rotor  stabilizer  bar. 
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VI.  CONCLUSION 


A,  MODEL  DISCREPANCIES 

The  current  model  can  only  perform  for  smaller  commanded  yaw  angles,  from 
0°to  about  20°.  Beyond  20°,  the  position  errors  get  larger  as  the  model  moves  through 
its  flight  path  and  will  even  enter  into  spiraling  motion  and  or  divergence  from  the 
desired  flight  path.  Further  tests  will  have  to  be  conducted  to  determine  the  cause  of  this 
unexpected  behavior.  The  tail  rotor  model  seems  to  be  doing  fine  on  its  own  given  that 
single  inputs  of  varying  yaw  angles,  with  yaw  controller,  generate  the  expected  yaw 
response.  A  step  in  yaw  angle  for  180°  and  270° both  gave  over  damped  yaw  angle 
responses  with  roughly  1 -second  settling  time.  Figure  53  illustrates  the  yaw  angle 
response  to  a  270°  step  command. 


Yaw  Step  vs  Angle  Response  (deg) 


Figure  53.  Yaw  Step  Response  t//'  =  270° 

Figure  54  illustrates  the  control  input  generated  by  the  PD  controller  in  order  to 
get  the  UAV  to  the  270°  orientation. 
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Figure  54.  Yaw  Step  Command  Control  Effort 

Figure  54  shows  that  the  eontrol  effort  is  only  required  while  getting  the  UAV  to 
the  desired  orientation,  from  1  to  roughly  2.5  seconds,  after  the  yaw  angle  has  settled. 
The  significance  of  the  control  signal  illustrated  in  Figure  54  is  that  it  shows  that  no 
control  effort  is  required  once  the  desired  orientation  is  obtained.  This  is  also  in  keeping 
with  the  fact  that,  once  the  desired  orientation  is  obtained,  and  control  effort  taken  out, 
the  tail  rotor  thrust  should  return  to  the  hover  value  of  7.29  Newtons,  as  discussed  in 
Chapter  III.C.2,  p.  42.  Figure  55  illustrates  the  tail  rotor  thrust  response  to  the  270°  step 
command. 


Control  Input  u 


1 - ^ - r 


J _ 1 _ . _ I _ I _ I _ 1 _ I _ L 


80 


Tail  Rotor  Thrust  Response  to  Step  Collective  (Newton) 
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Figure  55.  Tail  Rotor  Thrust  Response  to  Yaw  Step  Command/Statie 

At  1  seeond,  the  tail  rotor  thrust  is  driven  to  almost  -15  Newton  in  order  to 
generate  the  neeessary  yaw  moment  required  to  make  the  UAV  spin  to  +270° .  As  the 
vehiele  approaehes  the  desired  yaw  angle,  eounter-thrust  is  applied  to  keep  the  tail  from 
overshooting;  the  thrust  eontinues  to  oseillate  until  the  tail  rotor  thrust  eonverges  baek  to 
the  hovering  tail  rotor  thrust  of  7.29  Newton.  Figures  52  through  54  support  the  validity 
of  the  tail  rotor  model  design. 

The  problem  arises  when  the  yaw  angle  eommand  is  followed  up  with  a  position 
step  command  in  either  direction.  Figure  56  illustrates  the  position  response  to  the 
^ref  =  1  position  command.  The  position  step  command  is  given  4  seconds  after  the 

yaw  step  command  of  270° .  This  is  done  to  ensure  that  sufficient  settling  time  has  been 
given  to  the  yaw  angle. 
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Step  Position  Response  to  Xref=1  m 
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Figure  56.  Position  Response  at  =  270° 

By  the  time  the  model  reaches  the  270°  orientation,  the  coordinate  system  for  the 
model  should  be  oriented  such  that  the  inertial  frame  x  -axis  is  pointing  in  the  direction 
of  the  earth  frame’s  -y  -axis.  After  the  270°  turn,  the  step  command  in  the  x  direction 
should  result  in  the  model  flying  towards  the  body  fixed  x  direction  or,  the  earth 
coordinate  -y  direction.  What  Figure  56  illustrates  is  that  the  model  is  going  in  the 
direction  of  the  +x  earth  coordinate,  before  turning  back  into  the  -x  earth  coordinate 
direction. 

The  model  discrepancy  is  unlikely  to  be  due  to  tail  rotor  modeling.  Instead,  the 
most  likely  source  of  the  error  is  in  the  inertial  frame  direction  cosine  matrix  signals 
coming  from  the  6-DOF  blockset.  These  signals  will  have  to  be  analyzed,  and  signal 
gains  may  have  to  be  applied  in  order  to  obtain  the  expected  results. 
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All  of  the  simulations  using  the  box  and  Figure  9  are  sueh  that  the  earth 
eoordinates  and  inertial  frames  are  in  alignment  and  is  why  the  simulations  from  Chapter 
V  worked  fine  for  ^  =  0°  . 

B,  FUTURE  WORK 

Future  work  would  need  to  address  the  diserepaneies  between  the  body  fixed 
frame  and  the  earth  eoordinates,  assuming  that  this  is  the  souree  of  simulation  errors  for 
large  yaw  angle  eommands. 

Additional  future  work  ineludes  quantitative  model  verifieation  from  flight  test 
data  and  iterative  model  tuning  and  eonfiguration.  The  tuned  model  will  need  to  be  tested 
in  other  flight  eonfigurations  in  order  to  obtain  state  transition  matrix  derivative  funetions 
and  robust  eontrol  testing  will  need  to  be  applied.  Onee  the  model  has  been  fully 
validated  and  a  robust  eontroller  satisfaetorily  tested,  the  next  step  will  be  to  embed  the 
eontrol  teehnique  to  the  MieroPilot  Autopilot  for  Fiardware  in  the  Loop  testing. 
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APPENDIX  A:  MIT  INSTRUMENTED  X-CELL  60  SYSTEM 
IDENTIFICATION  UAV  PARAMETERS  (AFTER  [1]) 


m=8.2  kg 

helicopter  mass 

Ixx=.18  kg 

rolling  moment  of  inertia 

Iyy=.34  kg 

pitching  moment  of  inertia 

Izz=.28  kg 

yawing  moment  of  inertia 

Q=167  rad/sec 

Main  rotor  speed 

=5.5  rad  ‘ 

Main  rotor  blade  lift  curve  slope 

c;:=024 

Main  rotor  zero  lift  drag  coefficient 

Kr  =  -235  m 

Main  rotor  hub  height  above  c.g. 

l,^=.9\m 

Tail  rotor  hub  location  behind  eg 

=  .08  m 

Tail  rotor  height  above  eg 

=166  Nm/rad 

Longitudinal  flapping  spring  derivative  at  hover 

=  82.6  Nm/rad 

Longitudinal  flapping  spring  derivative  at  hover 
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APPENDIX  B:  NEWTON-RAPHSON  TECHNIQUE  FOR  ROTOR 
INFLOW  ESTIMATION  (M  FILES) 


getvi_v2,m 

%getvi_v2 . m 

%Determine  Main  Rotor  Inflow  using  Newton-Rhapson  Technique 
%Input:  in=[u  v  w  uwind  vwind  wwind  T  ro  alpha  vi[n+l]] 

%Output :  vi 

function  vi=getvi  v2(in) 
i=l; 

u=in  (1) ; 
v=in (2 ) ; 
w=in ( 3 ) ; 
uwind=in ( 4 )  ; 
vwind=in ( 5 ) ; 
wwind=in ( 6) ; 

T=in(7) ; 

ro=in ( 8 ) ; 

alpha=abs ( in ( 9 ) ) ; 

vi  previous=in ( 1 0 )  ; 

c=ro*  ((167*.775)''2)  *pi*  (.775''2)  ; 

CT=T/c; 

vt=sqrt  (  (u+uwind)  "^2+  (v+vwind)  "'2)  ; 
mu=vt / (167*. 775) ; 
vz=w+wwind; 
vh=4 . 171; 

lh=4.171/ (167*. 775)  ; 
mud=vz/ (167*. 775) ; 
if  vz<0  &&  vz>-.834 

lamda ( 1 ) =lh* ( 1-mud/ Ih) ; 
elseif  vz<=-.834  &&  vz>-2*4.171 

lamda(l)  =  (.11108*vz  +  5.0  97  6) / (167*.  7  7  5)  ; 
elseif  vz<=-2*4.171 

vi=-4.171* ( (vz/ (2*vh) ) +sqrt ( ( (vz/ (2*vh) ) ^2) -1) ) ; 
lamda(l)=vi/ (167*. 775)  ; 
elseif  vz>=0 

lamda (1) =4 . 171/ (167*.775) ; 

e=l; 

i=l; 

if  i<500 

while  abs (e) >=. 00005 

flamda (i) =lamda (i) -mu*tan (alpha) -CT/ (2*sqrt (mu^2+lamda (i) ^2) ) ; 

flamda  prime  (i)  =l+CT/2*  (  (mu"'2  +  lamda  (i)  ''2)  ''  (-3/2)  )  *lamda  (i)  ; 

lamda ( i+1 ) =lamda ( i ) -flamda ( i ) / flamda  prime ( i ) ; 

e=abs ( (lamda ( i+1 ) -lamda ( i ) ) / (lamda ( i+1 ) ) ) ; 

i=i+l ; 

end 

elseif  i>=500 

lamda (i) =vi_previous/ (167*  .  775)  ; 
end 

end 

vi=lamda(i) *167*. 775; 
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getvi_v2TR,m 


%getvi_v2TR.m 

%Determine  Tail  Rotor  Inflow  using  Newton-Rhapson  Technique 
%Input:  in=[u  v  w  uwind  vwind  wwind  Ttr  ro  alpha=0  vi[n+l]] 

%Output :  vi 

function  vi=getvi  v2(in) 
i=l; 

u=in  (1) ; 
v=in (2 ) ; 
w=in ( 3 ) ; 
uwind=in ( 4 )  ; 
vwind=in ( 5 ) ; 
wwind=in ( 6) ; 

T=in(7) ; 

ro=in ( 8 ) ; 

alpha=abs ( in ( 9 ) ) ; 

c=ro* ( (77  8. 2*. 13) -2) *pi*  (.13^2)  ; 

CT=T/c; 

vt=sqrt  (  (u+uwind)  "^2+  (v+vwind)  "'2)  ; 
mu=vt/ (778. 2*. 13)  ; 
vz=w+wwind; 
vh=7 .294; 

lh=7 .294/ (778.2* . 13) ; 
mud=vz/ (778.2*. 13)  ; 
if  vz<0  &&  vz>-.834 

lamda ( 1 ) =lh* ( 1-mud/ Ih)  ; 
elseif  vz<=-.834  &&  vz>-2*4.171 

lamda (l)  =  (.11108*vz  +  5. 0976) / (778. 2*. 13); 
elseif  vz<=-2*4.171 

vi=-7.294* ( (vz/ (2*vh) ) +sqrt ( ( (vz/ (2*vh) ) ^2 ) -1 ) ) ; 
lamda (1) =vi/ (778. 2*. 13); 
elseif  vz>=0 

lamda (1) =7 .294/ (778. 2*. 13); 

e=l; 

i=l; 

if  i<500 

while  e> . 0005 

flamda  (i)  =lamda  (1)  -mu*tan  (alpha)  -CT/  (2*sqrt  (mu''2  +  lamda  (i)  ''2)  ) 
flamda  prime  (i)  =l+CT/2*  (  (mu"'2  +  lamda  (i)  ''2)  ''  (-3/2)  )  *lamda  (i)  ; 
lamda ( i+1 ) =lamda ( i ) -flamda ( i ) / flamda  prime ( i ) ; 
e=abs (lamda ( i+1 ) -lamda ( i ) ) / (lamda ( 1+1 ) ) ; 

1=1+1 ; 
end 

elseif  i>=500 

lamda (1) =7 .294/ (778. 2*. 13); 
end 

end 

vi=lamda (1) *778.2* . 13; 
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UpdateMR,m 


%UpdateMR 

%  1  2  3  4  5  6  7  8  9  10  11  12  13 

%Inputs:u=[u  v  w  uwind  vwind  wwind  ro  vt  cp  cp[n+l]  T[n+1]  vi 
vi [n+1] ] 

function  out=updateMR ( in) 

u=in  (1) ; 

v=in (2 ) ; 

w=in ( 3 ) ; 

uwind=in ( 4 )  ; 

vwind=in ( 5 ) ; 

wwind=in ( 6) ; 

ua=u+uwind; 

va=v+vwind; 

ro=in ( 7 ) ; 
vt=in ( 8 ) ; 
cp=in (9) ; 

cp_previous=in  (10)  ; 

T  previous=in ( 1 1 ) ; 
vi=in ( 12 ) ; 
li=vi/ (167*. 775) ; 
vi  previous=in (13) ; 

li  previous=vi  previous/ (167* . 775) ; 

fx= . 1 ; 
fy=.22; 

Dx=fx*ro*ua*vt/ 2 ; 

Dy=f y*ro*va*vt/ 2 ; 

D=sqrt ( (Dx) ^2+ (Dy) ^2) ; 
alpha=(57.3*D*pi/180) /  (80.414) ; 

%Define  Input  to  thrust  function 

%input_f_thrust= [u  v  w  uwind  vwind  wwind  vi  cp  ro  cp[n+l]  T[n+1]] 
input  f  thrust=[u  v  w  uwind  vwind  wwind  vi  cp  ro  cp  previous 
T_previous]  ; 

%output:  Thrust 

T=getThrustMR ( input  f  thrust); 

%Define  Input  to  get  induced  velcity  function 
%Input:  in=[u  v  w  uwind  vwind  wwind  T  ro  alpha  vi[n+l]] 
input  f  vi=[u  V  w  uwind  vwind  wwind  T  ro  alpha  vi  previous]; 
vi=getvi  v2 (input  f  vi); 
out= [T  vi] ; 
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APPENDIX  C:  PARAMETER  DEFINITION  FILES 


Config.m 


%Fixed  Parameters 


mass=8 . 2 ; 
g=  9 . 8  0  6  6  5 ; 

Inertia=[.18  0  0;0  .34  0;0  0  .28]  ; 
ICParam. Ts=l / 1 00 ;  %Sample  time  period 
SR=1 / ICParam. Ts ;  %Sample  Rate 

g _ 

o - 

%Flapping  Dynamics  Constants 

0 _ 

O - 


LN=3 . 92; 
omega=l 67 ; 

KB=18 . 89729; 

IB=. 038; 
tauf= . 1 ; 

g _ 

o - 

%Main  Rotor  Moment 

g _ 

o - 


%Lock;  Number 
%MR  speed  rad/sec 

%hub  torsional  stiffness  {NM/rad} 
%mr  blade  flapping  inertia  {kg  m''2} 


Derivatives 


Ma=82 . 6; 


Lb=166; 


%Vehicle  Geometric  Properties 


hmr= .235; 
Svf=.012; 
Rtr=. 13; 
Ltr= . 91 ; 
htr=. 08 


%height  of  MR  from  eg  position 
%Vertical  fin  area 
%tail  rotor  radius 
%tail  hub  loc  behind  eg 
%tr  height  above  eg 


9- 

o 


%Transmitter  Duty  Cycle  Limits 
%[.6E-3  1.52E-3  2.4E-3]  sec  duty  cycle 
%[-l  0  1]  input  in  stick 

g _ 

o - 


corresponds 


to 


duty  cycle  limits= [ . 6E-3  1.52E-3  2.4E-3]; 
lat_stick_positions= [-1  0  1]; 
lon_stick_positions= [-1  0  1]; 
pedal_positions= [-1  0  1]; 
collective_stick_positions= [ 0  .5  1.0]; 
duty_cycle_col= [ . 6E-3  .9E-3  2.4E-3]; 
servo_output= [-45  0  45]; 

MR  collective  range= [-3 . 88*pi/180  5.495*pi/180  15*pi/180]; 

TR  collective  range= [-20 . 0*pi/180  1 1 . 22 91 *pi/ 1 8 0  25*pi/180]; 
max  cyclic  command=20; 

%  Name  of  the  MAT-file  that  will  be  generated 
cfgmatfile  =  ' TestPlantcfg ' ; 

%  Save  workspace  variables  to  MAT  file 
save (cfgmatfile) ; 

%  Output  a  message  to  the  screen 

fprintf ( streat ( ' \n  Aircraft  configuration  saved  as:\t', 
streat (cfgmatfile) , ' .mat ' ) ) ; 
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fprintf ( ' \n' ) ; 


CreateModelStructure.m 

%CreateModel structure .  m 

0 _ 

O - 

%Default  Initial  Conditions:  Hovering 

g _ 

o - 

ICParam.pos= [0  0  0]; 

ICParam.vel= [0  0  0]; 

ICParam. euler= [ 0  0  0]; 

ICParam.pqr= [0  0  0]; 

ICParam. flap  angle=[0  0]; 

ICParam. MRThrust=80 . 414; 

ICParam. MRInflow=4 .171; 

ICParam. MRCollectivePitch=5 . 495*pi/180; 

ICParam. TRInf low=7 .294; 

ICParam. TRThrust=6 . 92  ; 

ICParam. YawMoment=0 . 0; 

ICParam. ro=l .225; 

ICParam. Wind= [0  0  0]; 

%Simulink;  Model  to  Trim 

%ICParam. SimModel=input (' Enter  name  of  model  to  run  and  trim:') 

ICParam. SimModel= ' NonLinearTRexAlign ' ; 

fprintf (' \nThe  model  order  is  being  determined ... \n ') ; 

%Get  the  sim  options  structure 

%To  run  simget,  ensure  the  model  is  set  to  fixed  step  in  the 
%Simulation>>Conf iguration  Parameters  dialogue,  with  the  step  time  set 
%to  ICParam. Ts  defined  above 

runtime=input (' Enter  the  Simulation  Run  time  (sec)\n'); 

ICParam. SimOptions=simget (ICParam. SimModel) ; 

LatCyc=0 ; 

LonCyc=0 ; 

TRCol=ll .2291; 

MRCol=5 .495; 

Trimlnput= [LonCyc  LatCyc  MRCol  TRCol] *pi/180; 

%Run  Model  for  1  sample  period 

fprintf (' Please  wait,  running  simulation  at  hovering  conditions...') 
[SimTime, SimStates , SimOutputs ] =sim ( ICParam. SimModel , [0  . 01 ] , . . . 

ICParam. SimOptions, [0  Triminput; . 01  Triminput] ) ; 

%Find  the  state  order 

InitialStates= [ ICParam. pos ' ; ICParam. euler ' ; ICParam. vel ' ; ICParam. pqr ' ; . . 
ICParam. flap  angle']; 

ICParam. NHeliStates=length ( InitialStates ) ; 

ICParam. NSimulinkStates=length (SimStates (1,  : ) )  ; 

ICParam. StateIdx=zeros (ICParam. NHeliStates, 1) ; 
clear  SimTime  SimStates  SimOutputs 
fprintf ( ' Done . \n ' )  ; 

clear  SimTime  SimStates  SimOutputs; 
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APPENDIX  D:  TRIMMING  PROCEDURE  FILE 


TrimHeli.m 


%TrimHeli . m 

%%%Determine  The  Initial  Guess  For  Aircraft  Controls 
%Choose  Flight  Condition  to  trim 

fprintf (' Enter  1  to  trim  at  Hover,  2  to  trim  at  other  specific 
statesXn ' ) 

condition=input (' Enter  your  choice : \n ') ; 

if  condition==l 

ICParam.vel= [0  0  0]; 

ICParam.pos= [0  0  0]; 

ICParam. euler= [ 0  0  0]; 

ICParam. Wind= [0  0  0]; 
elseif  condition==2 
fprintf ( ' \n' ) ; 

fprintf (' \nChoose  Steady  State  Flight  Condition  to  Trim:'); 

fprintf  (  '  \n - \n  '  ) 

ICParam. vel=input (' Enter  Desired  Velocity  States  [u  v  w] : ' ) ; 

ICParam. pos=input (' Enter  Desired  Position  [x  y  z]:'); 

ICParam. euler=input (' Enter  Desired  Euler  Angles  [fi  theta  psi]  in 
deg: ' ) *pi/180; 

ICParam. Wind=input (' Enter  the  Wind  Vector  to  apply  to  the  model  [uw  vw 

ww]  : ' )  ; 

%Note,  pqr  and  flap  angles  must  be  zero  in  steady  state 
end 

%Trim  Error  Threshold 

MaxErrVel= [ . 05  .05  .05];  %Not  more  than  5cm/sec  in  error 
MaxErrEuler= [ 1  1  l]*pi/180;  %Not  more  than  .25deg  in  angle  error 
%MaxErrVel= [ . 2  .2  .2]; 

%MaxErrEuler= [ 3  3  3]*pi/180; 

%Control  Surface  Gains 
vel=max (ICParam. vel) ; 

K=-.5/15*vel+l; 

GainLat=- . 01 ; 

GainLon= . 01 ; 

GainTR=. 01; 

GainCollective=- . 01 ; 

fprintf ( ' \n  Computing  the  initial  estimates  for  the  required  trim 
inputs . . . \n ' ) ; 

GoodGuess=0;  Niter=l; 

Gain= [GainLon  GainLat  GainCollective  GainTR] ; 

goodguess=[0  0  0  0]; 

while  (GoodGuess==0) && (Niter<50) 

%Velocity  Loop 

%Run  Simulink  Model  for  a  short  time  (10  sec) 

fprintf (' Please  wait,  running  simulation  and  comparing  simulation 
states  with  desired  states...') 
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[ SimTime, SimStates , SimOutputs ] =sim ( I C Par am. SimModel , [0  2] . 

ICParam. SimOptions,  [0  Triminput;  2  Triminput] ) ; 
ErrVel=SimOutputs (end, 4 : 6) -ICParam. vel; 

ErrEuler=SimOutputs (end, 7:9) -ICParam. euler; 
fprintf ( ' \nIteration  #%2d\n' , Niter) 
fprintf  (  ' - \n'  ) 

fprintf (' \nX  Velocity  Error=%2.3f  m/s\n ' , ErrVel ( 1 ) ) 

fprintf('Y  Velocity  Error=%2.3f  m/s\n ', ErrVel (2 ) ) 

fprintf (' Z  Velocity  Error=%2.3f  m/s\n ', ErrVel (3) ) 

fprintf (' \nRoll  Angle  Error=%2.3f  (deg) \n ' , ErrEuler ( 1 ) *  1 8 0/pi ) 

fprintf ( '  Pitch  Angle  Error=%2.3f  (deg) \n ', ErrEuler (2 )* 1 8 0/pi ) 

fprintf ( '  Yaw  Angle  Error=%2.3f  (deg) \n ', ErrEuler (3) *180/pi) 

for  i=l:3; 

if  abs (ErrVel ( i )) <=MaxErrVel ( i ) 

%TrimInput (i) =TrimInput (i) ; 
goodguess ( i ) =1 ; 

else 

Triminput (i) =TrimInput (i) +ErrVel (i) *Gain (i)  ; 

end 

end 

if  abs (ErrEuler (3) ) <=MaxErrEuler (3) 

%TrimInput ( 4 ) =TrimInput ( 4 ) ; 
goodguess ( 4 ) =1 ; 

else 

Triminput ( 4 ) =TrimInput ( 4 ) +ErrEuler ( 3 ) *Gain ( 4 ) ; 

end 

Niter=Niter+l ; 

if  goodguess== [ 1  111] 

GoodGuess=l ; 

else 

GoodGuess=0 ; 

exit=input  (  '  Hit  1  to  Escape  Loop,  Enter  to  continue:'),' 
if  exit==l 

GoodGuess=l ; 

else 

GoodGuess=0 ; 

end 

end 

end 

%Save  Initial  Guesss 

%?Trim  condition  for  Ion  and  lat  for  constant  velocity? 
TrimParam.MRCol=TrimInput (3) ; 

TrimParam.TRCol=TrimInput (4) ; 

TrimParam. Lon=TrimInput ( 1 )  ; 

TrimParam. Lat=TrimInput (2 ) ; 

TrimParam. vel=SimStates (end,  7:9)  '  ; 

TrimParam. pos=SimStates (end,  1:3)  '  ; 

TrimParam. euler=SimStates (end,  4:6)  ' ; 

TrimParam. f laps=SimStates (end,  13:14)  '  ; 

TrimParam. pqr=SimStates (end,  10:12)  '  ; 
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clear  SiinTime  SimStates  SimOutputs  Triminput 


fprintf ( ' \nlnitial  Guesses  for  trim  Inputs  are:  Servo  Linkage 
Command (deg) \n ' ) 

fprintf  (  ' -  - 

- \n '  ) 

fprintf ( '  MR  Collective  Command 

%2 . 3f \n ' , TrimParam.MRCol*180/pi) 
fprintf ( '  TR  Collective  Command 

%2 . 3f \n '  ,  TrimParam. TRCol*18  0/pi) 
fprintf (' Longitudinal  Flap  Angle  Command 
%2 . 3f \n ' , TrimParam. Lon* 18  0/pi) 
fprintf ( '  Lateral  Flap  Angle  Command 
%2 . 3f \n ' , TrimParam. Lat* 180/pi) 

%Perform  Helicopter  Trim 
fprintf ('\n'); 

fprintf (' \nPerforming  the  aircraft  trim...\n') 

%Set  Initial  Guesses 

InitialStates=zeros (ICParam.NSimulinkStates, 1) ; 

InitialStates (1) =TrimParam.pos (1)  ; 

InitialStates (2) =TrimParam.pos (2)  ; 

InitialStates (3) =TrimParam.pos (3)  ; 

InitialStates (4) =TrimParam. euler ( 1 )  ; 

InitialStates (5) =TrimParam. euler (2) ; 

InitialStates (6) =TrimParam. euler (3)  ; 

InitialStates (7) =TrimParam. vel ( 1 )  ; 

InitialStates (8) =TrimParam. vel (2) ; 

InitialStates (9) =TrimParam. vel (3) ; 

InitialStates (10) =TrimParam.pqr ( 1 ) ; 

InitialStates (11) =TrimParam.pqr (2) ; 

InitialStates (12) =TrimParam.pqr (3)  ; 

InitialStates (13) =TrimParam. flaps (1)  ; 

InitialStates ( 14 ) =TrimParam. flaps (2 ) ; 

%  Set  optimization  parameters 

TrimParam. Options ( 1 )  =  1;  %  show  some  output 

TrimParam. Options ( 14 )  =  1000;  %  max  iterations 

Initiallnput= [TrimParam. Lon; TrimParam. Lat; TrimParam. MRCol; TrimParam. TRC 
ol]  ; 

InitialOutput= [TrimParam. pos ; TrimParam. vel ; TrimParam. euler] 
InitialDerivatives=zeros (ICParam.NSimulinkStates,  1)  ; 

%StateFixIdx  is  a  list  of  indeces  for  the  states  that  need  to  be 
%held  constant.  In  this  case,  [phi  theta  psi  u  v  w  alpha  beta  p  q  r] 
StateFixIdx=[4  5  6  7  8  9  10  11  12  13  14]; 

InputFixIdx= [ ] ; 

OutputFixIdx= [ 4  56789]; 

%State  derivatives  to  be  held  fixed 
DerivFixIdx=[7  8  9  10  11  12  13  14]; 

%Trim  The  Helicopter 

[TrimOutput . States , TrimOutput . Inputs , TrimOutput . Outputs . . . 
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, TrimOutput . Derivatives , options ] =trim ( I C Par am. SimModel , Initial States 
, Initialinput, InitialOutput, StateFixIdx, InputFixIdx, OutputFixIdx 
, InitialDerivatives , DerivFixIdx, TrimParam. Options) ; 
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APPENDIX  E:  RUN  MODEL  FILE 


%RunModel . m 

%This  program  allows  the  user  to  choose  desired  states,  run  the  model 
with 

%the  chosen  states,  find  the  trim  conditions  for  the  nonlinear  model, 
find 

%the  linear  representation,  and  plot  step  responses  from  a  central 

command 

%program 

clear  all 

clf 

%Run  the  configuration  file  that  defines  all  the  model  parameters 
run  Config 

run  CreateModelStructure 

exit=0 ; 

trimed=0 ; 

linear=0 ; 

while  exit==0 


fprintf ( ' \nMENU\n ' ) 

fprintf('You  must  Choose  Option  1  Before  any  Other  Choice\n') 
fprintf ( '================================================\n' ) 


fprintf ( 
fprintf ( 
fprintf ( 
fprintf ( 
fprintf ( 
fprintf ( 
fprintf ( 
fprintf ( 
fprintf ( 


1- 

2- 

3- 

4- 

5- 

6- 

7- 

8- 
9- 


Choose  Desired  States  &  Trim  Heli  Model\n') 

Print  Trim  Results\n') 

Generate  Step  Responses  for  NonLinear  Model\n') 
Plot  Step  Responses  for  Linear  Model  &  Compare\n') 
Display  Heave  Model  Transfer  Functions\n ' ) 


Display  the  Yaw  Dynamics  Model 
Display  the  Lon  and  Lat  Linear 
Conduct  Derivative  Analysis\n') 
Exit\n ' ) 


Transfer  Function\n') 
Models  Derived\n') 


menu=input (' Enter  Your  Choice :\n') 


if  menu==l 

run  TrimHeli 
trimed=l ; 

fprintf (' NonLinear  Model  Has  Been  Trimmed\n') 
elseif  menu==2 

if  trimed==0 

fprintf ('You  must  First  Choose  States  and  Trim,  Choose  1  at 

the  MenuXn ' ) 

elseif  trimed==l 

%Print  the  Trim  Results 

fprintf ( ' \n=================================\n ' ) ; 

fprintf (' \nThe  Trim  Results  Are:'); 

fprintf ( ' \n - 


(deg) 

fprintf ( ' \n 
' , TrimOutput . Inputs ( 3 ) 

MR 

*180/pi) ; 

Collective 

Command  = 

%3.3f 

(deg) 

fprintf ( ' \n 
' , TrimOutput . Inputs  ( 4 ) 

TR 

*180/pi) ; 

Collective 

Command  = 

%3.3f 

(deg) 

fprintf ( ' \n 
' , TrimOutput . Inputs ( 1 ) 

Longitudinal 
*180/pi) ; 

Cyclic  Command  = 

%3.3f 
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fprintf ( ' \n  Lateral  Cyclic  Command  =  %3.3f 

(deg) \n ' , TrimOutput . Inputs (2)  *180/pi)  ; 
else 

exit=0 ; 

end 

elseif  menu==3 
if  trimed==0 

fprintf ('You  must  First  Choose  States  and  Trim,  Choose  1  at 

the  MenuXn ' ) 

elseif  trimed==l 

clear  UT  Triminput 
time=0: (1/SR) :runtime; 
samples=runtime*SR+l ; 

TrimInput=zeros (4, samples)  ; 

Linear . SimModel= ' LinearTRexAlign ' ; 

Linear. SimOptions=simget (Linear. SimModel) ; 

for  i=l:4 

%TrimInput ( 1 , : ) =TrimParam. Lon; 

%TrimInput (2 , : ) =TrimParam. Lat; 

%TrimInput (3, : ) =TrimParam.MRCol; 

%TrimInput ( 4 , : ) =TrimParam. TRCol ; 

Triminput ( i , : ) =TrimOutput . Inputs ( i ) ; 

end 

size=input (' Enter  desired  step  size  in  degrees : \n ') ; 

Step Input (1, 101 : samples) =pi/180*size; 

LonStep=TrimInput ( 1 , : ) ' fStepInput ' ; 

LatStep=TrimInput (2 , : ) ' IStepInput ' ; 

MRColStep=TrimInput ( 3 ,  :)  ' fStepInput ' ; 

MRColStepDown=TrimInput ( 3 ,  :)  ' -Stepinput ' ; 

TRCol Step=TrimInput ( 4 , ; ) ' + Step Input ' ; 

UTl=[time'  LonStep  Triminput (2,  :)  '  Triminput (3,  :)  ' 

Triminput (4,  : )  '  ]  ; 

UT2=[time'  Triminput  (1,  :)  '  LatStep  Triminput  (3,  :)  ' 

Triminput (4,  : )  '  ]  ; 

UT3up=[time'  Triminput (1,  :)  '  Triminput (2,  :)  '  MRColStep 
Triminput  (4,  :  )  '  ]  ; 

UT3dwn= [time '  Triminput (1,  :)  '  Triminput (2,  :)  '  MRColStepDown 
Triminput (4,  : )  '  ]  ; 

UT4=[time'  Triminput  ( 1 ,:)  '  Triminput  (2 ,;)  '  Triminput  ( 3 ,:)  ' 

TRColStep] ; 

fprintf (' Please  wait,  running  simulation  ... \n ' ) 

[ SimTime,  SimStates ,  SimOutputsl ] =sim ( I C Par am. SimModel ,  [ 0 

runtime] , . . . 

ICParam. SimOptions, UTl) ; 
t=SimTime; 

figure  ( 1 ) 

ymin=min (LonStep) -.25*pi/180; 
ymax=max (LonStep) +.25*pi/180; 
subplot ( 3 , 1 , 1 ) 

plot  (t, UTl ( : , 2) *180/pi) ,  grid  on,  axis([0  runtime 
ymin*180/pi  ymax*180/pi] ) 

title (' Longitudinal  Step  Input  (deg) ') 
subplot ( 3 , 1,2) 
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plot ( t , SimOutputsl ( : , 4 ) , ' bo ' , t , SimOutputsl ( : , 2 ) , t , SimOutputsl ( : , 3 ) , ' g-- 

'  ) 

grid  on,  title (' Velocity  (m/s) ') 

legend ('u', 'v', 'w', ' location ' , ' northwest ' ) ; 

subplot (3,1,3) 

plot  (t, SimOutputsl (:,8)*180/pi,  'r- 
' , t, SimOutputsl (:,8)*180/pi,  'bo') 

grid  on,  title ('Euler  Angles  (deg) ') 

legend ( ' theta ' , ' theta  linear ' , ' location ' , ' northwest ' ) ; 
fprintf ( ' Done  Running  Longitudinal  Step  Simulation\n ' ) 

[ SimTime, SimStates , Sim0utputs2 ] =sim ( I C Par am. SimModel , [ 0 

runtime] , . . . 

ICParam. SimOptions, UT2) ; 
figure (2 ) 

ymin= (min (LatStep) -.25*pi/180) *180/pi; 
ymax= (max (LatStep) +.25*pi/180) *180/pi; 
subplot  ( 3 , 1 , 1 ) 

plot (t, UT2 ( : , 3) *180/pi) ,  grid  on,  axis([0  runtime  ymin 

ymax] ) 

title (' Lateral  Step  Input  (deg) ') 
subplot ( 3 , 1,2) 

plot ( t ,  Sim0utputs2 ( : , 4 ) , t , Sim0utputs2  ( : , 5 ) ,  ' bo ' , t , Sim0utputs2 ( : , 6) ,  ' g-- 

'  ) 

grid  on,  title (' Velocity  (m/s) ') 

legend ('u',  'v',  'w',  ' location ' ,  ' northwest ' ) ; 

subplot (3,1,3) 

plot (t, Sim0utputs2 (:,7)*180/pi,  'r- 
' , t, Sim0utputs2 (:,8)*180/pi,t, Sim0utputs2 (;,9)*180/pi,  'bo') 
grid  on,  title ('Euler  Angles  (deg) ') 
legend ( ' phi ' , ' theta ' , ' psi ' , ' location ' , ' northwest ' ) ; 

[SimTime, SimStates , Sim0utputs3up] =sim ( ICParam. SimModel , [ 0 

runtime] , . . . 

ICParam. SimOptions, UT3up) ; 

%Extract  heave  dynamics  climbing  model  from  data  obtained 
from  nonlinear 

%model 

cpMR=MRColStep; 
w=Sim0utputs3up ( : , 6) ; 
heavedata=iddata (w, cpMR, 1/SR) ; 
heavemodel  climb=n4sid (heavedata, 3) ; 

[e, xoMRc] =pe (heavemodel  climb, heavedata) ; 
fprintf ('The  Climb  Heave  Model  is:\n') 

[ahc, bhc, chc, dhc] =ssdata (heavemodel_climb) ; 

[numh,  denh] =ss2tf (ahc,  bhc,  chc,  dhc) ; 
hc=tf (numh, denh) ; 

figure  ( 3 ) 

ymin= (min (MRColStep) -.25*pi/180) *180/pi; 
ymax= (max (MRColStep) +.25*pi/180) *180/pi; 
subplot (2 , 1 , 1 ) 
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ymax] ) 

plot (t, UT3up ( : , 4) *180/pi)  ,  grid  on,  axis([0  runtime  ymin 

title ('MR  Step  Input  (deg) ') 
subplot  (2 , 1,2) 

plot (t, SimOutputs3up ( : , 6) , ' bo ' ) 
grid  on,  title (' Velocity  (m/s) ') 
legend ( ' w ' , ' location ' , ' northwest ' ) ; 

runtime] , . 

[ SimTime, SimStates , SimOutputs3dwn] =sim ( I C Par am. SimModel , [ 0 

ICParam. SimOptions, UT3dwn) ; 

%Extract  heave  dynamics  descent  model  from  data  obtained 
from  nonlinear 

%model 


cpMRd=MRColStepDown; 
wd=SimOutputs3dwn ( : , 6) ; 
heavedatad=iddata (wd, cpMRd, 1/SR) ; 
heavemodel  descent=n4sid (heavedatad, 3) ; 

[e, xoMRd] =pe (heavemodel  descent, heavedatad) ; 
fprintf('The  Descent  Heave  Model  is:\n') 

[ahd, bhd, chd, dhd] =ssdata (heavemodel  descent) ; 

[numd, dend] =ss2tf (ahd, bhd, chd, dhd) ; 
hd=tf (numd, dend) ; 

ymax] ) 

figure ( 4 ) 

ymin= (min (MRColStepDown) -.25*pi/180) *180/pi; 
ymax= (max (MRColStepDown) +.25*pi/180) *180/pi; 
subplot (2 , 1 , 1 ) 

plot  (t, UT3dwn ( : , 4) *180/pi) ,  grid  on,  axis([0  runtime  ymin 

title ('MR  Step  Input  (deg) ') 
subplot (2 , 1,2) 

plot (t, SimOutputs3dwn ( : , 6)  ,  'bo ' ) 
grid  on,  title (' Velocity  (m/s) ') 
legend ( ' w ' , ' location ' , ' northwest ' ) ; 

runtime] , . 

[SimTime, SimStates , SimOutputs4 ] =sim ( ICParam. SimModel , [ 0 

ICParam. SimOptions, UT4) ; 

model 

cpt=UT4 ( : , 5 ) ;  %tail  rotor  input  data 
psi=SimOutputs4 ( : , 9) ; 

yawdata=iddata (psi, cpt, 1/SR) ;  %create  the  IDDATA  object 
yawmodel=n4sid (yawdata) ;  %generate  the  identification 

[e, xot] =pe (yawmodel, yawdata) ;  %obtain  error  and  initial 

states 

fprintf('The  Yaw  Dynamics  TF  is:\n') 

[ay, by, cy, dy] =ssdata (yawmodel) ; 

[numy, deny] =ss2tf (ay,  by,  cy,  dy)  ; 
yaw=tf (numy, deny) ; 

figure ( 5 ) 
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ymin= (min (TRColStep) -.25*pi/180) *180/pi; 
ymax= (max (TRColStep) +.25*pi/180) *180/pi; 
subplot (2 , 1 , 1 ) 

plot  (t, UT4 ( : , 5) *180/pi)  ,  grid  on,  axis([0  runtime  ymin 

ymax] ) 

title ('TR  Step  Input  (deg) ') 
subplot (2 , 1,2) 

plot (t,  Sim0utputs4 (:,9)*180/pi,  'bo') 
grid  on,  title ('Euler  Angles  (deg) ') 
legend ( ' psi ' , ' location ' , ' northwest ' ) ; 


end 

elseif  menu==4 

fprintf  (  '  Note  that  the  model  that  will  be  extracted  uses  the 
current  trim  settings . \n ' ) 

%Extract  the  linear  model 

fprintf (' \n  \nExtracting  Helicopter  Linear  Model ... \n ') ; 
%Perturbation  Level 
LinParam  (1)  =10''-8; 

[A, B, C, D] =linmod ( ICParam. SimModel , TrimOutput . States ,  .  .  . 

TrimOutput . Inputs , LinParam) ; 

ssplant=ss (A, B, C, D) 

run  ExtractSISO 

linear=l ; 

fprintf ('Lat  and  Long  Model  Extracted  at  chosen  operation 
pointXn ' ) 

[ SimTime, LinStates , LinOutputsl ] =sim (Linear . SimModel ,  [ 0 
runtime] , . . . 

ICParam. SimOptions, UTl) ; 
figure  ( 1 ) 

ymin= (min (LonStep) -.25*pi/180) *180/pi; 
ymax= (max (LonStep) +.25*pi/180) *180/pi; 
subplot ( 3 , 1,2) 

plot ( t , SimOutputsl ( : , 4 ) , ' bo ' , t , SimOutputsl ( : , 5 ) , t , SimOutputsl ( : , 6) , ' g-- 
'  , t , LinOutputsl ( : , 4 ) ,  ' g .  ' ) 

legend ('u',  'v',  'w',  'u  linear ' ,  ' location ' ,  ' northwest ' ) ; 
grid  on 

title (' Velocity  (m/s)') 
subplot (3,1,3) 

plot ( t , SimOutputsl (:,8)*180/pi, 'bo',t, LinOutputsl (:,8)*180/pi, 'g.') 
grid  on,  title ('Euler  Angles  (deg) ') 

legend ( ' theta ' , ' theta  linear ' , ' location ' , ' northwest ' ) ; 

[SimTime, LinStates , Lin0utputs2 ] =sim (Linear . SimModel ,  [ 0 
runtime] , . . . 

ICParam. SimOptions, UT2) ; 

figure (2 ) 
subplot ( 4 , 1,2) 

plot (t, Sim0utputs2 ( : , 4 ) , ' b . ' , t, Sim0utputs2 ( : , 5) , 'bo ' , t, Sim0utputs2 ( : , 6) 

,  ' r-- ' , t, Lin0utputs2 ( :  ,  5)  ,  ' g.  '  ) 
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grid  on,  title (' Velocity  (m/s) ') 

legend ('u', 'v', 'w', 'v  linear ' , ' location ' , ' northwest ' ) ; 
subplot (4,1,3) 

plot ( t , Sim0utputs2 (:,7)*180/pi, 'bo',t, Lin0utputs2 (:,7)*180/pi, 'g.') 
grid  on,  title ('Euler  Angles  (deg) ') 
legend ( ' phi ' , ' phi  linear ' , ' location ' , ' northwest ' ) ; 
subplot (4,1,4) 

plot (t, Sim0utputs2 ( : , 9) *180/pi) 

[ SimTime, LinStates , Lin0utputs3up] =sim (Linear . SimModel , [ 0 
runtime] , . . . 

ICParam. SimOptions, UT3up) ; 

%Extract  heave  dynamics  climbing  model  from  data  obtained  from 
nonlinear 

%model 

cpMR=UT3up ( : , 4) ; 
heavedata=iddata (w, cpMR, 1/SR) ; 
heavemodel  climb=n4sid (heavedata,  3) ; 

[e, xoMRc] =pe (heavemodel  climb, heavedata) ; 
fprintf('The  Climb  Heave  Model  is:\n') 

[ahc, bhc, chc, dhc] =ssdata (heavemodel_climb) ; 

[numh, denh] =ss2tf (ahc,  bhc,  chc,  dhc) ; 
hc=tf (numh, denh) 

figure ( 3 ) 
subplot (2 , 1,2) 

plot (t, Sim0utputs3up ( : , 6) , ' bo ' , t, Lin0utputs3up ( : , 6) , ' g . ' ) 

grid  on,  title (' Velocity  (m/s) ') 

legend ( ' w '  ,  ' w  linear '  ,  ' location '  ,  ' northwest ' ) ; 

[SimTime, LinStates , Lin0utputs3dwn] =sim (Linear . SimModel ,  [ 0 
runtime] ,... Trimlnput= [LonCyc  LatCyc  MRCol  TRCol] *TrimInput= [LonCyc 
LatCyc  MRCol  TRCol] *pi/180;pi/180; 

ICParam. SimOptions, UT3dwn)  ; 

cpMR=UT3dwn ( :  ,  4 )  ; 

heavedata=iddata ( Sim0utputs3dwn ( : ,  6)  ,  cpMR,  1/SR)  ; 
heavemodel_descent=n4sid (heavedata,  3)  ; 

[e, xoMRd] =pe (heavemodel_descent,  heavedata)  ; 
fprintf('The  Descent  Heave  Model  is:\n') 

[ahd, bhd, chd, dhd] =ssdata (heavemodel_descent) ; 

[numd, dend] =ss2tf (ahd, bhd,  chd,  dhd)  ; 
hd=tf (numd, dend) 

figure ( 4 ) 
subplot (2 , 1,2) 

plot ( t , Sim0utputs3dwn ( ; ,  6 )  ,  ' bo ' , t , Lin0utputs3dwn ( : ,  6)  ,  ' g . '  ) 

grid  on,  title (' Velocity  (m/s) ') 

legend ( ' w ' , ' w  linear ' , ' location ' , ' southwest ' ) ; 

[SimTime, LinStates , Lin0utputs4 ] =sim (Linear . SimModel , [ 0 
runtime] , . . . 
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ICParam. SimOptions, UT4) ; 


%extract  linear  model  for  yaw  dynamics  using  SysID  toolbox 
%Extract  the  Yaw  dynamics  model  from  data  obtained  from 
nonlinear 

%model 

cpt=UT4 ( : , 5 ) ;  %tail  rotor  input  data 

yawdata=iddata (psi, cpt, 1/SR) ;  %create  the  IDDATA  object 
yawmodel=n4sid (yawdata) ;  %generate  the  identification 

model 

[e, xot] =pe (yawmodel, yawdata) ;  %obtain  error  and  initial  states 
fprintf('The  Yaw  Dynamics  TF  is:\n') 

[ay, by, cy, dy] =ssdata (yawmodel) ; 

[numy, deny] =ss2tf (ay, by, cy,  dy) ; 
yaw=tf (numy, deny) 

figure ( 5 ) 
subplot (2 , 1,2) 

plot ( t , Sim0utputs4 (:,9)*180/pi,  'bo',t, Lin0utputs4 (:,9)*180/pi,  'g.') 
grid  on,  title ('Euler  Angles  (deg) ') 
legend ( 'psi' , 'psi  linear' , 'location' , ' southwest ' ) ; 
elseif  menu==5 

fprintf('The  Climbing  Heave  Model  is:\n') 
he 

fprintf('The  Descent  Heave  Model  is:\n') 
hd 

elseif  menu==6 

fprintf('The  Yaw  Dynamics  Model  is:\n') 
yaw 

elseif  menu==7 

fprintf('The  Longitudinal  Models  are:\n') 
fprintf ( ' dlon  to  x\n') 

Ion  X  tf 

fprintf (' dlon  to  u\n') 

Ion  u  tf 

fprintf (' dlon  to  theta\n') 

Ion  theta  tf 
pause 

fprintf ('The  Lateral  Models  are:\n') 

fprintf (' dlat  to  y\n') 

lat_y_tf 

fprintf (' dlat  to  v\n') 
lat_v_tf 

fprintf (' dlat  to  phi\n') 
lat  phi  tf 
elseif  menu==8 

%Resize  Matrix  A  to  fit  Padfield  pg.  210 

%go  from  state  x=[x  y  z  phi  theta  psi  u  v  w  p  q  r  alpha  beta] 
to 

%  x=[u  w  q  theta  v  p  phi  r] 
dummyA=zeros (8) ; 
dummyB=zeros (8,4) ; 
column=[7  9  11  5  8  10  4  12]; 
count=0 ; 
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for  j  =1 : 8 
count=count+l ; 
for  i=l : 8 

col=column ( i ) ; 

dummyA ( j , i) =A (column ( j ) 

if  count==l 

dummyB (i,  : ) =B (col,  :  ) ; 

else 

end 


end 


col ) 


end 

Apadf ield=dummyA 
Bpadf ield=dummyB 

elseif  menu==9 
exit=l ; 

end 

end 
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APPENDIX  F:  EXTRACT  SINGLE  INPUT  SINGLE  OUTPUT 

MODELS 


%Extract  SISO  Models 

%Mak;e  sure  you  already  have  the  linear  plant  model  defined 
%The  model  ssplant  is  defined  in  RunModel.m,  menu  option  3 

%Define  the  linear  longitudinal  model:  input=dlon  output=[x 
%Model  Inputs=[dlon  dlat  MRCol  TRCol] 

%  1234567  8  9 

%ModelOutputs= [x  y  z  u  v  w  phi  theta  psi] 


0 _ 

O - 

%Longitudinal  Model 

0 _ 

O - 

Ion  to  x=minreal ( ssplant ( 1 ,  1 ))  ; 
lon_to_y=minreal (ssplant (2,1)); 

Ion  to  z=minreal ( ssplant ( 3 ,  1 ))  ; 

Ion  to  u=minreal ( ssplant ( 4  ,  1 ))  ; 

Ion  to  v=minreal ( ssplant ( 5 ,  1 ))  ; 

Ion  to  w=minreal (ssplant (6,  1)  )  ; 

Ion  to  phi=minreal ( ssplant ( 7  ,  1 ))  ; 
lon_to_theta=minreal (ssplant (8,1) ) ; 
lon_to_psi=minreal (ssplant (9,1)  )  ; 


g _ 

o - 

%Lateral  Model 

g _ 

o - 

lat  to  x=minreal ( ssplant ( 1 , 2  ))  ; 
lat_to_y=minreal (ssplant (2,2)  )  ; 
lat_to_z=minreal (ssplant (3,2)  )  ; 
lat_to_u=minreal (ssplant (4,2)  )  ; 
lat_to_v=minreal (ssplant (5,2)  )  ; 
lat_to_w=minreal (ssplant (6,2)  )  ; 
lat  to  phi=minreal ( ssplant ( 7 , 2  )  )  ; 
lat_to_theta=minreal (ssplant (8,2) ) ; 
lat_to_psi=minreal (ssplant (9,2)  )  ; 


g _ 

o - 

%Heave  Dynamics  Model 

g _ 

o - 

MRCol  to  x=minreal ( ssplant ( 1 , 3 ))  ; 
MRCol_to_y=minreal (ssplant (2,3))  ; 
%MRCol_to_z=minreal (ssplant (3,3)  )  ; 
MRCol  to  u=minreal ( ssplant ( 4 , 3 ) )  ; 
MRCol  to  v=minreal ( ssplant ( 5 , 3 )  )  ; 
%MRCol_to_w=minreal (ssplant (6,3)  )  ; 
MRCol  to  phi=minreal ( ssplant ( 7 , 3 ) )  ; 
MRCol_to_theta=minreal (ssplant (8,3)  )  ; 
MRCol_to_psi=minreal (ssplant (9,3) )  ; 


u  theta] 
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%Yaw  Dynamics  Model 

0 _ 

O - 

TRCol  to  x=minreal ( ssplant ( 1 , 4 )  )  ; 
TRCol_to_y=minreal (ssplant (2,4)); 
TRCol  to  z=minreal ( ssplant ( 3 , 4 ) ) ; 
TRCol  to  u=minreal ( ssplant ( 4 , 4 ) )  ; 
TRCol  to  v=minreal ( ssplant ( 5 , 4  )  )  ; 
TRCol  to  w=minreal (ssplant (6,  4)  )  ; 
TRCol  to  phi=minreal ( ssplant ( 7 , 4 ) )  ; 
TRCol_to_theta=minreal (ssplant (8,4)  )  ; 
%TRCol_to_psi=minreal (ssplant (9,4) ) ; 
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APPENDIX  H:  MINIMUM  REALIZATION  TRANSFER  FUNCTION 
MODELS  AT  HOVER  (1°  STEP  COMMAND) 


Longitudinal  Dynamics  Transfer  Functions 

X  2.842x10 ‘"5' +6.821x10 +98.055^+33695' +2.405x10"5  +  1633  ,  , 

- = - ^ ^ - 2 - ^ ^ -  (m/rad) 

5® +44.395' +590.35" +24885' +243.55' +5. 1895-3.628x10^*' 


u  _  2.842x10 ‘"5" +98.055' +33695' +2.405x10"5  +  1633 
5'  +44.395"  +590.35'  +24885'  +243.55  +  5.189 


(m/sec-rad) 


6>  3.553x10 ‘"5" +7.958x10 ‘'5' -0.064285' -0.64285  +  1.231x10^’  , 

- = - ^ - 2 - ^ ^ -  (rad/rad) 

5' +44.295" +585.95' +24295' +0.00032495 +  0.001347 


Lateral  Dynamics  Transfer  Functions 


y  _  -5.4x10 ‘'5'  -7.071x10"5"  +98.055'  +1.003x10"5'  +9.074x10"5  +  2817 
5*’ +112.35' +19565" +94I65' +923.65' +19.75  +  6.944x10  *' 

(m/rad) 


V 


Plat 


-1.421x10  *'5"  +98.055'  +1.003x10"5'  +9.074x10"5  +  2817 
5'  +112.35"  +19565'  +94I65'  +923.65  +  19.7 


(m/sec-rad) 


A. 

Plat 


2.842x10  ‘"5"  -4.547x10  ‘'5'  +0.12145'  +1.2145  +  1.671x10^^ 
5'  +112.25" +  19445'  +92225'  +0.0010785  +  0.005114 


(rad/rad) 


Climbing  Heave  Dynamics  Transfer  Functions 

w  2.O6I5' -2.735 +  .6697  ,  , 

—  =  ^ ^ - r  (m/sec -rad) 

6>^  5' -1.9855' +.98545-5.7x10  " 


Descending  Heave  Dynamics  Transfer  Functions 


w 

P 


2.0645' -3.9835-1.919 
5' -2.93 15' +2.8625-0.9317 


(m/sec-rad) 
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Yaw  Dynamics  Transfer  Function 


^  _  1.691x10^^5^ -0.0128/ +0.01615^0.0089535^ -0.01615  +  0.00383 

5^-3.6685' +4.625^-1.715' -0.93665^0.8445-0.1545  ) 


MODE 

EIGENVALUES 

Longitudinal  Velocity  (u) 

[-24.3  -10.14  -9.86  -.0686  -.0312] 

Pitch 

[-24.3  -10.3  -9.7  ±7xl0V] 

Lateral  Velocity  (v) 

[-92.19  -10.3  -9.7  -.0685  -.0312] 

Roll 

[-92.2  -10±.1567  ±7.44x10^V] 

no 


LIST  OF  REFERENCES 


[1]  Vladislav  Gavrilets,  “Autonomous  Aerobatic  Maneuvering  of  Miniature 
Helicoptersrotary-wing,"  PhD  Thesis,  MIT,  Boston,  MA,  2003. 

[2]  Bernard  Mettler,  Identification  Modeling  and  Characteristics  of  Miniature 
Rotorcraft,  Kluwer  Academic  Publishers,  2003. 

[3]  Raymond  W.  Prouty,  Helicopter  Performance,  Stability,  and  Control,  PWS 
Publishers,  1986. 

[4]  Gareth  D.  Padfield,  Helicopter  Flight  Dynamics:  The  Theory  and  Application  of 
Flying  Qualities  and  Simulation  Modeling,  AIAA  Education  Series,  1996. 

[5]  J.  Gordon  Leishman,  Principles  of  Helicopter  Aerodynamics,  Cambridge 
University  Press,  2000. 

[6]  John  D.  Anderson,  Jr.,  Fundamentals  of  Aerodynamics,  McGraw-Hill,  Inc.,  1991. 

[7]  Jim  Ledin,  Embedded  Control  Systems  in  C/C++,  CMP  Books,  2004. 

[8]  Katsuhiko  Ogata,  Modern  Control  Engineering,  Pearson  Prentice  Hall,  2003. 

[9]  Katsuhiko  Ogata,  MATLAB  for  Control  Engineers,  Pearson  Prentice  Hall,  2008. 

[10]  Norman  S.  Nise,  Control  Systems  Engineering,  John  Wiley  &  Sons,  Inc.,  2004. 


Ill 


THIS  PAGE  INTENTIONALLY  LEET  BLANK 


II2 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Teehnieal  Information  Center 
Ft.  Belvoir,  VA 

2.  Dudley  Knox  Library 
Naval  Postgraduate  Sehool 
Monterey,  CA 

3.  Dr.  Jeffrey  Knorr 

Eleetrieal  and  Computer  Engineering  Department 
Code  EC/Ko 

Naval  Postgraduate  Sehool 
Monterey,  CA 

4.  Dr.  Robert  G.  Hutchins 

Electrical  and  Computer  Engineering  Department 
Naval  Postgraduate  School 
Monterey,  CA 

5.  Dr.  Vladimir  Dobrokhodov 
Mechanical  Engineering  Department 
Naval  Postgraduate  School 
Monterey,  CA 

6.  loannis  Kitsios 

Mechanical  Engineering  Department 
Naval  Postgraduate  School 
Monterey,  CA 

7.  Marine  Corps  Representative 
Naval  Postgraduate  School 
Monterey,  CA 

8.  Director,  Training  and  Education,  MCCDC,  Code  C46 
Quantico,  VA 

9.  Director,  Marine  Corps  Research  Center,  MCCDC,  Code  C40RC 
Quantico,  VA 

10.  Marine  Corps  Tactical  Systems  Support  Activity  (Attn;  Operations  Officer) 
Camp  Pendleton,  CA 


113 


